Numerical Investigation of Elastic Layer Effects in Wheel–Rail Rolling Contact

: In railway systems, layered structures could be induced in wheel–rail contact interfaces due to several causes, such as head hardening, work hardening, plastic deformation, and mechanical or thermal excursion-induced phase transformation. This study proposes an explicit ﬁnite element (FE) method for investigating elastic layer effects in wheel–rail rolling contact. The proposed method is ﬁrst validated by comparing its solution with that of Kalker’s boundary element method (BEM) when the layer is not present, with a focus on the tractive rolling contact. To investigate general layer effects, the rail is assumed to consist of two layers, i.e., the top layer and the matrix material. The top layer is assumed to have different elastic moduli from the matrix material and then the top elastic layer effects on contact characteristics such as contact stress, contact patch, and subsurface stress are investigated. Different layer thicknesses are also considered. It is observed that a harder layer tends to introduce larger contact pressure and surface shear stress, but a smaller contact patch. A harder layer also produces larger subsurface stresses. A thicker layer may intensify these effects. The results suggest that in engineering applications, the analysis of wheel–rail rolling contact consequences such as wear and rolling contact fatigue (RCF) may need to consider the layered structures using appropriate methods.


Introduction
The wheel-rail rolling contact problem is a crucial issue in railway systems.The calculation of rolling contact characteristics is important for solving many contact problems, including wear, plastic deformation, RCF, vehicle-track interaction, and traction control.In many previous studies, the analysis of these contact problems used traditional halfspace-based methods such as Vermeulen and Johnson's analytical method [1] and Kalker's BEM [2], which assume that the contact bodies are half-spaced and have isotropic elasticity.However, in reality, the material properties are complicated and the structures could be layered.Therefore, to accurately evaluate if it is necessary to consider the layer in further applications, it is crucial to study the layer effects on contact characteristics, including stresses at the interface and stresses in the subsurface.
In railway systems, layered structures present in the wheel-rail contact interfaces are unavoidable.For instance, head-hardened rails are widely installed in railway networks to increase wear resistance and fatigue strength [3].The tensile strength and hardness of headhardened rails can be significantly increased compared with standard carbon rails, from 880 MPa and 260 HB to 1130 MPa and 340 HB for head layer material, respectively [4].In addition, while the rail is in service, the rail surface is affected by mechanical and thermal excursions under repeated wheel-rail rolling contact loads.As a result, the topmost surface layer of the rail is subjected to cyclic plastic deformation and phase transformation, where the microstructure evolves and its property changes to form a harder layer or a white etching layer (WEL) [5].The plastic deformation layer can reach a thickness of 1 mm.A WEL is brittle and has a hardness of over 700 HV [3].The thickness of the WEL can be up to 2 mm [6].Meanwhile, the elastic modulus increases with the hardness [7].Moreover, when repairing the rail with laser cladding, a cladding material is added to the rail surface, producing a bonded clad layer on the rail surface [8].This means that a layered structure is formed on the rail's top surface.Also, third bodies in railway engineering can be considered as layered structures [9].In these cases, the rails become layered structures or gradient structures.Besides rails, layered structures are also observed on the surfaces of wheels.For instance, through the metallurgical observation of wheel tread damages in the field, Zhang et al. [10] found that a microstructure change could occur in the wheel surface and become large amounts of upper bainite.This upper bainite had higher hardness and elasticity than the matrix microstructure.This layered structure was associated with stress concentration and RCF crack initiation.Zhou et al. [11] also found that inhomogeneities in materials can affect rolling contact fatigue life in rotating mechanical components.Additionally, Du et al. [12] conducted experimental and numerical studies and found that multi-layered materials can induce tensile failure in rolling contact systems, indicating that layered structures exist in rails and wheels.How they influence the wheel-rail rolling contact should be evaluated.
Some effort has been made to solve rolling contact problems associated with layers.Since the 1960s, analytical and numerical models have been developed to study the effect of elastic and viscoelastic layers bonded to an elastic half-space in rolling contact.For instance, Burton [13] proposed an analytical method for solving visco-elastic effects in the lubrication of rolling contact.Margetson [14] analyzed the rolling contact of a rigid cylinder over a smooth elastic or viscoelastic layer.The bonded and non-bonded layer effects on contact pressure were also studied.Londhe et al. [15] proposed an extended Hertz theory for the contact of case-hardened steels.In these above studies, only normal contact problems were considered, and tangential problems were not included.To solve the tangential problems associated with layers, Kalker [16] proposed a numerical method based on the BEM for two-dimensional visco-elastic multilayered cylinders rolling contact.To study the influence of a third body layer, Meierhofer et al. [17] proposed a model for calculating the traction characteristic for the wheel-rail contact with a third body layer.However, in these studies, only two-dimensional problems were considered.Analytical and numerical methods were also developed to solve three-dimensional rolling contact problems with layers.Zhang and Yan [18] proposed a semi-analytical model for the effects of subsurface inhomogeneous on frictional rolling contact.Xi et al. [19] recently proposed a numerical model for solving three-dimensional rolling contact problems with elastic coating layers.Goryacheva and Miftakhova [20] developed a model for an elastic sphere rolling over a thin viscoelastic layer bonded to an elastic half-space.Guler et al. [21] used a numerical method for the tractive rolling contact problem between a rigid cylinder and a graded coating.Although some of these studies considered both normal and tangential problems associated with layers, they assumed the substrate to be half-spaces.This means that the layers are bonded or not bonded to the half-spaces.The contact geometries were also simplified.
In reality, more complicated material properties and actual geometries have to be considered in the wheel-rail contact as they have significant effects on wheel-rail contact behavior; thus, the half-space assumption has to be abandoned.FE methods can meet the requirements.To this end, an explicit FE method has been developed to solve rolling contact problems.The Explicit FE method is suitable for solving wheel-rail transient rolling contact [22].With this method, all possible creepages, complex materials, realistic geometries, and dynamic effects can be included.It was first applied to the dynamic wheel-rail frictional rolling contact by Zhao et al. [22].Later, it was used for rolling contact with spin [23] and compression-shift-rolling contact with partial slip [24].Its accuracy has also been validated against traditional methods by approximating it to quasi-static states.Furthermore, this method has been applied to more complex geometries and materials such as surface defects [25], corrugation [26], and material discontinuities [27].In the dynamic wheel-rail contact, it was compared with methods based on multibody dynamics [28] and was shown to have better performance.All these studies focused on the wheel-rail rolling contact and wheel-rail dynamic interactions with isotropic materials.However, the effects on the wheel-rail rolling contact characteristics have not been particularly studied using the explicit FE method where the layered structures of the rail are involved.
In this study, the explicit FE method proposed in [22] is employed to investigate the influence of rail surface layered structure on wheel-rail rolling contact characteristics such as contact stresses, contact patch, and stick-slip.The study considers a general case of the gradient layered structure instead of focusing on a particular engineering problem such as the head-hardened rail surface or WEL.The top layer is assumed to have a different elastic modulus from the matrix material.Different elastic moduli and layer thicknesses are considered.In total, five thicknesses, from 0.2 mm to 6 mm, are included; four different elastic moduli are employed for each thickness.The resulting stick-slip condition, surface contact stresses, and subsurface stresses in the rail are investigated.

Problem Statement
The wheel-rail rolling contact problem with a layered rail structure is schematically described in Figure 1.Each layer is assumed to be homogeneous and linearly elastic, the thickness of the layer is h, and the elastic moduli for layer h and the matrix material are represented by E 1 and E 0 , respectively.The corresponding Poisson's ratios are v 1 and v 0 .The wheel material is assumed to be the same as the base rail material.In this case, the layer effect on the wheel-rail rolling contact characteristics is studied.Semi-analytical/analytical models are not applicable to this case due to the complicated geometries of the wheel and rail.
Lubricants 2023, 11, x FOR PEER REVIEW 3 of 16 [23] and compression-shift-rolling contact with partial slip [24].Its accuracy has also been validated against traditional methods by approximating it to quasi-static states.Furthermore, this method has been applied to more complex geometries and materials such as surface defects [25], corrugation [26], and material discontinuities [27].In the dynamic wheel-rail contact, it was compared with methods based on multibody dynamics [28] and was shown to have better performance.All these studies focused on the wheel-rail rolling contact and wheel-rail dynamic interactions with isotropic materials.However, the effects on the wheel-rail rolling contact characteristics have not been particularly studied using the explicit FE method where the layered structures of the rail are involved.In this study, the explicit FE method proposed in [22] is employed to investigate the influence of rail surface layered structure on wheel-rail rolling contact characteristics such as contact stresses, contact patch, and stick-slip.The study considers a general case of the gradient layered structure instead of focusing on a particular engineering problem such as the head-hardened rail surface or WEL.The top layer is assumed to have a different elastic modulus from the matrix material.Different elastic moduli and layer thicknesses are considered.In total, five thicknesses, from 0.2 mm to 6 mm, are included; four different elastic moduli are employed for each thickness.The resulting stick-slip condition, surface contact stresses, and subsurface stresses in the rail are investigated.

Problem Statement
The wheel-rail rolling contact problem with a layered rail structure is schematically described in Figure 1.Each layer is assumed to be homogeneous and linearly elastic, the thickness of the layer is h, and the elastic moduli for layer h and the matrix material are represented by  and  , respectively.The corresponding Poisson's ratios are  and  .The wheel material is assumed to be the same as the base rail material.In this case, the layer effect on the wheel-rail rolling contact characteristics is studied.Semi-analytical/analytical models are not applicable to this case due to the complicated geometries of the wheel and rail.

FE Model of Wheel-Rail Rolling Contact
An FE model of wheel-rail rolling contact is constructed to study the layer effects.A rail and a typical wheel with primary suspension and mass block of a car body are included in this model, as illustrated in Figure 2. Since the wheel-rail contact in a steady state is investigatedfixed, the track and vehicle structures are not included in the model.The related parameters are from [29].An 8000 kg mass is applied.The gravity of the mass is transmitted to the wheel through primary suspension.The primary suspension is

FE Model of Wheel-Rail Rolling Contact
An FE model of wheel-rail rolling contact is constructed to study the layer effects.A rail and a typical wheel with primary suspension and mass block of a car body are included in this model, as illustrated in Figure 2. Since the wheel-rail contact in a steady state is investigatedfixed, the track and vehicle structures are not included in the model.The related parameters are from [29].An 8000 kg mass is applied.The gravity of the mass is transmitted to the wheel through primary suspension.The primary suspension is applied to minimize the vibration between the wheel and rail.Its stiffness and damping are 0.95 MN/m and 10.4 kN s/m, respectively.The weight of the wheel is 660 kg.LMA and CN60 are employed as the wheel and rail profile types, respectively.The wheel-rail rolling contact problem can be considered a 3D transient structural dynamic problem.The dynamic problem is widely described by constrained partial differential governing equations [30], which can be transformed into semi-discrete equilibrium equations for the solution using an explicit integration scheme with a central difference.This solution procedure is implemented in ANSYS/LS-DYNA.A surface-to-surface contact searching scheme is applied in this method and a penalty method is employed to enforce the kinematic normal contact condition [31].The explicit FE method also includes the computation of tangential contact, which is dealt with using the Coulomb friction law.
The FE model is meshed with 3D 8-node solid elements.To ensure precision and to reduce the computation cost at the same time, a non-uniform mesh is employed in the model.Therefore, a fine mesh of 1 mm is assigned to the potential contact surfaces and the investigated area in subsurface, as shown in Figure 2. A coarser mesh is applied to other regions.The mass block is modeled using eight mass elements and the primary suspension is modeled using eight groups of spring and damping elements.The parameters of the mass and suspension are averagely distributed into the eight groups of elements in the model.
When rolling the wheel, a tangential load is often produced between the wheel and rail during driving or braking.In this study, the driving condition is considered.In this process, the tangential force is produced through the transmission of the torque applied to the motor.In the model, a torque,  , is directly applied on the wheel axis.Consequently, a tangential force,  , is produced between the wheel and the rail.
In the simulation, the wheel is set to roll over the rail from left to right along the x direction.The simulation consists of two steps [17].The first step involves an implicit analysis to calculate the initial static equilibrium state of the system under gravity.The second step involves explicit analysis with stress initialization.The results from the first step are initially imported into the explicit analysis to initialize the system.The explicit analysis starts at time zero with an initialized system and an applied initial velocity.A torque is applied on the wheel axle to drive the wheel.In this case, the wheel-rail contact is a The wheel-rail rolling contact problem can be considered a 3D transient structural dynamic problem.The dynamic problem is widely described by constrained partial differential governing equations [30], which can be transformed into semi-discrete equilibrium equations for the solution using an explicit integration scheme with a central difference.This solution procedure is implemented in ANSYS/LS-DYNA.A surface-to-surface contact searching scheme is applied in this method and a penalty method is employed to enforce the kinematic normal contact condition [31].The explicit FE method also includes the computation of tangential contact, which is dealt with using the Coulomb friction law.
The FE model is meshed with 3D 8-node solid elements.To ensure precision and to reduce the computation cost at the same time, a non-uniform mesh is employed in the model.Therefore, a fine mesh of 1 mm is assigned to the potential contact surfaces and the investigated area in subsurface, as shown in Figure 2. A coarser mesh is applied to other regions.The mass block is modeled using eight mass elements and the primary suspension is modeled using eight groups of spring and damping elements.The parameters of the mass and suspension are averagely distributed into the eight groups of elements in the model.
When rolling the wheel, a tangential load is often produced between the wheel and rail during driving or braking.In this study, the driving condition is considered.In this process, the tangential force is produced through the transmission of the torque applied to the motor.In the model, a torque, M T , is directly applied on the wheel axis.Consequently, a tangential force, F T , is produced between the wheel and the rail.
In the simulation, the wheel is set to roll over the rail from left to right along the x direction.The simulation consists of two steps [17].The first step involves an implicit analysis to calculate the initial static equilibrium state of the system under gravity.The second step involves explicit analysis with stress initialization.The results from the first step are initially imported into the explicit analysis to initialize the system.The explicit analysis starts at time zero with an initialized system and an applied initial velocity.A torque is applied on the wheel axle to drive the wheel.In this case, the wheel-rail contact is a tractive rolling contact.The tractive effort is determined by a traction coefficient, µ, which is defined by the static gravity expressed as where F T is the tangential force exerted by the torque M T and G is the static weight of the wheel and car body.In this study, two traction coefficients, 0.15 and 0.2, are employed in the validation with the friction coefficient of 0.5.The 0.2 traction coefficient is used to study the layer effects, associated with the same friction coefficient.The traction and friction coefficients can be varied in future studies.The rolling contact solution is obtained when the wheel reaches position 0.3 on the x-axis, indicating that the wheel has reached a steady state.Thus, a quasi-static solution is obtained.

Simulation Cases
In this study, both the wheel and the rail are assumed to be elastic.The layer effects of the rail are investigated.The layered surface is produced by setting Young's modulus to be different from the matrix material.The layer and the matrix material are bonded.A Young's modulus of 210 GPa is applied to E 0 .E 0 and a Poisson's ratio of 0.3 are used for the matrix material of the wheel and rail.A density of 7800 kg/m 3 is applied to both the surface layer and the matrix material.In the case without the layer, E 1 has the same value as E 0 .To study the effect of the layer, three elastic moduli (E 1 = 1.67E 0 , E 1 = 2E 0 and E 1 = 4E 0 ) are compared with those for the case without the layer.Other parameters are the same as those of the matrix material.Different layer thicknesses (0.2, 0.5, 1, 3, 6 mm) are modeled.The width of the layer is set to 60 mm in the lateral direction and the contact patch in this direction is located roughly in the middle of the width.

Process of Contact Solutions
The contact pressure, surface shear stress, stick and slip areas, and subsurface stresses produced during wheel-rail contact are investigated in this study.The subsurface stresses are directly obtained from the post-process software LS-PrePost of ANSYS/LS-DYNA.The contact pressure, surface shear stress, and stick and slip areas are determined by processing the calculated nodal forces as in [23].
A node is in the contact patch if where F n_N is the nodal force in the direction normal to the contact surface, and is calculated from the nodal force in three axis directions: F x , F y , and F z .The surface tangential force is the total projection of the nodal forces in three-axis directions.In the current study, the frictional force is approximately equal to the nodal force in the z-axis direction.The contact pressure and the surface shear stress are obtained by dividing the normal nodal force F n_N and the frictional force F f with the area of the face on the rail surface of the block element, respectively.The stick and slip areas are distinguished by comparing the limiting frictional force with the frictional force.If the frictional force on a node satisfies then the node is in the stick area.Otherwise, the node is in the slip area.

Validation
Although the explicit FE method for rolling contact problems has been validated in [22], the mesh size in the present study is different from those used in previous studies.To determine whether the mesh used in the present study is reasonable, this section validates the accuracy of the FE model by comparing the solutions with those of Kalker's BEM [2] when the layer is not included, as shown in Figure 3.In the comparison, the contact pressure and surface shear stress along the central line of the contact patch of the tractive rolling contact with two traction coefficients are investigated.The differences in the maximum stresses measured using the two methods are within 2%; thus, both the contact pressure and the surface shear stresses from the two methods are in good agreement.However, unlike the maximum stress, a difference of nearly 18.7% is found in the size of the contact patch.This difference may arise due to different treatments in the two methods, as discussed in [22].A finer mesh could reduce the difference; however, the accuracy in the present study is usually acceptable in engineering applications.Additionally, the mesh size should not influence the layer effects in the present study since it is constant in all simulated cases.
pressure and surface shear stress along the central line of the contact patch of the tractive rolling contact with two traction coefficients are investigated.The differences in the maximum stresses measured using the two methods are within 2%; thus, both the contact pressure and the surface shear stresses from the two methods are in good agreement.However, unlike the maximum stress, a difference of nearly 18.7% is found in the size of the contact patch.This difference may arise due to different treatments in the two methods, as discussed in [22].A finer mesh could reduce the difference; however, the accuracy in the present study is usually acceptable in engineering applications.Additionally, the mesh size should not influence the layer effects in the present study since it is constant in all simulated cases.

Layer Effects on Contact Stresses at the Interface
This section investigates the effect of the surface layer elastic modulus on rolling contact behavior when a top layer is present.Three surface layer elastic moduli ( = 1.67 ,  = 2 and  = 4 ) are compared with the case without a layer.Different layer thicknesses (0.2, 0.5, 1, 3, and 6 mm) are also considered.The 0.2 traction coefficient is used.The solutions are analyzed when the wheel center is at a position 0.3 m on the x-axis.

Layer Effect on Contact Pressure
Figure 4 shows the calculated pressure as a function of both the elastic modulus and the thickness of the top layer.As a fact, the pressure distribution is three-dimensional.To easily observe the results, the pressure along the central line of the contact patch in the travel direction is presented.Four elastic moduli are investigated for each thickness.When  =  at a thickness of 0.2 mm (Figure 4a), there is no layer, and the wheel and rail materials are homogenous.In this case, the rolling contact is approximately equivalent to the Hertzian type.The solution was validated against the Hertz theory and Kalker's BEM [22,23].The contact pressure distribution is nearly parabolic, with a maximum pressure of 1075 MPa.When the elastic modulus of the top layer is changed, the rail becomes a layered structure.As the elastic modulus is increased, the layer becomes harder.The contact pressure distribution still exhibits a parabolic shape, whereas the maximum pressure increases slightly with the increase in elastic modulus.An increase of 1.4% is observed when the elastic modulus is increased from  to 4 .
When the thickness of the layer is increased, the effect of the elastic modulus on the contact pressure becomes more significant; the contact pressure becomes larger and the contact patch becomes shorter with an increase in the elastic modulus, as shown in Figure

Layer Effects on Contact Stresses at the Interface
This section investigates the effect of the surface layer elastic modulus on rolling contact behavior when a top layer is present.Three surface layer elastic moduli (E 1 = 1.67E 0 , E 1 = 2E 0 and E 1 = 4E 0 ) are compared with the case without a layer.Different layer thicknesses (0.2, 0.5, 1, 3, and 6 mm) are also considered.The 0.2 traction coefficient is used.The solutions are analyzed when the wheel center is at a position 0.3 m on the x-axis.

Layer Effect on Contact Pressure
Figure 4 shows the calculated pressure as a function of both the elastic modulus and the thickness of the top layer.As a fact, the pressure distribution is three-dimensional.To easily observe the results, the pressure along the central line of the contact patch in the travel direction is presented.Four elastic moduli are investigated for each thickness.When E 1 = E 0 at a thickness of 0.2 mm (Figure 4a), there is no layer, and the wheel and rail materials are homogenous.In this case, the rolling contact is approximately equivalent to the Hertzian type.The solution was validated against the Hertz theory and Kalker's BEM [22,23].The contact pressure distribution is nearly parabolic, with a maximum pressure of 1075 MPa.When the elastic modulus of the top layer is changed, the rail becomes a layered structure.As the elastic modulus is increased, the layer becomes harder.The contact pressure distribution still exhibits a parabolic shape, whereas the maximum pressure increases slightly with the increase in elastic modulus.An increase of 1.4% is observed when the elastic modulus is increased from E 0 to 4E 0 .
ever, the effect is not linear.For instance, in the case of the 3 mm thickness, the change in the maximum contact pressure is not apparent when the elastic modulus is increased from 2 0 to 4 0 .
In summary, the elastic modulus of the top surface layer affects the contact pressure.In general, there is a tendency for harder layers to induce larger pressures.However, the degree of the effect depends on the thickness of the layer.

Layer Effect on Surface Shear Stress
Tractive rolling contact induces tangential contact behavior, resulting in surface shear stress.The produced surface shear stresses corresponding to the aforementioned cases are given in Figure 5.The stress distribution is shown along the central line of the contact patch.Surface shear stress is also a function of the elastic modulus and the thickness of the top layer.As shown in Figure 5a, the surface shear stress for the case without the layer is a typical partial slip type.In the slip region (left side of point A), the distribution shape follows that of the pressure because the shear stress is the friction coefficient times of the pressure based on Coulomb's law.In the stick region (right side of point A), the distribution suddenly drops.In the case of the 0.2 mm thickness, the material becomes a layered structure as the elastic modulus increases, the surface shear stress exhibits a similar distribution shape, and the maximum value is increased accordingly.The maximum pressure increases to nearly 5% when  1 is increased from  0 to 4 0 .
Figure 5b-e shows that the layer effects become more significant as the thickness is increased.A harder layer tends to increase the maximum shear stress in all cases, although fluctuations are observed.When the thickness is increased from 0.2 to 6 mm, the increase in the maximum surface shear stress is nearly 12%.This indicates that a thicker layer When the thickness of the layer is increased, the effect of the elastic modulus on the contact pressure becomes more significant; the contact pressure becomes larger and the contact patch becomes shorter with an increase in the elastic modulus, as shown in Figure 4d,e.Compared with the 0.2 mm thickness, an increase of 14.4% in the maximum pressure is observed in the case of 6 mm thickness.This indicates that the effect of the elastic modulus on the contact pressure is also associated with the thickness of the surface layer.However, the effect is not linear.For instance, in the case of the 3 mm thickness, the change in the maximum contact pressure is not apparent when the elastic modulus is increased from 2E 0 to 4E 0 .
In summary, the elastic modulus of the top surface layer affects the contact pressure.In general, there is a tendency for harder layers to induce larger pressures.However, the degree of the effect depends on the thickness of the layer.

Layer Effect on Surface Shear Stress
Tractive rolling contact induces tangential contact behavior, resulting in surface shear stress.The produced surface shear stresses corresponding to the aforementioned cases are given in Figure 5.The stress distribution is shown along the central line of the contact patch.Surface shear stress is also a function of the elastic modulus and the thickness of the top layer.As shown in Figure 5a, the surface shear stress for the case without the layer is a typical partial slip type.In the slip region (left side of point A), the distribution shape follows that of the pressure because the shear stress is the friction coefficient times of the pressure based on Coulomb's law.In the stick region (right side of point A), the distribution suddenly drops.In the case of the 0.2 mm thickness, the material becomes a layered structure as the elastic modulus increases, the surface shear stress exhibits a similar distribution shape, and the maximum value is increased accordingly.The maximum pressure increases to nearly 5% when E 1 is increased from E 0 to 4E 0 .
induces a larger surface shear stress.The influence on the surface shear stress could be similar to that on pressure because the surface shear stress is proportional to the pressure in the slip area based on Coulomb's law.In conclusion, the shear stress tends to increase with the increase in the elastic modulus, but the degree of this effect also depends on the thickness of the layer and the value of the elastic modulus.

Layer Effects on Stick-Slip and Size of Contact Patch
The stick-slip condition is one of the most important phenomena in tractive rolling contacts.Therefore, this section investigates the layer effects on the stick-slip condition.The effect on the size of the contact patch is also analyzed.
Figure 6 shows the stick-slip distributions of the contact patch for different elastic moduli when the layer thickness is 0.2 mm.The green nodes are in the slip area, and the blue nodes are in the stick area.It can also be observed that the size of the contact patch tends to increase as the elastic modulus increases, although very slightly.Moreover, the shapes of the slip and stick regions are very similar, but the ratio increases from 56.6% to 59.1% when the elastic modulus changes from  0 to 4 0 .This means that the ratio of the slip area tends to increase as the layer becomes harder.
To compare the stick-slip conditions and the contact patch sizes for the case with the 0.2 mm thickness, Figure 7 gives the stick-slip distribution of the contact patch for the four different elastic moduli when the layer thickness is 6 mm.The ratio of the slip area gradually reduces from 56.6% for  0 to 50.0% for 4 0 .At the same time, the contact patch becomes smaller.This could be because the harder layer produces smaller deformation when the layer becomes thicker, resulting in a smaller slip area and a smaller contact patch.In this case, the contact pressure and surface shear stress become larger, as shown in Figures 4 and 5  Figure 5b-e shows that the layer effects become more significant as the thickness is increased.A harder layer tends to increase the maximum shear stress in all cases, although fluctuations are observed.When the thickness is increased from 0.2 to 6 mm, the increase in the maximum surface shear stress is nearly 12%.This indicates that a thicker layer induces a larger surface shear stress.The influence on the surface shear stress could be similar to that on pressure because the surface shear stress is proportional to the pressure in the slip area based on Coulomb's law.In conclusion, the shear stress tends to increase with the increase in the elastic modulus, but the degree of this effect also depends on the thickness of the layer and the value of the elastic modulus.

Layer Effects on Stick-Slip and Size of Contact Patch
The stick-slip condition is one of the most important phenomena in tractive rolling contacts.Therefore, this section investigates the layer effects on the stick-slip condition.The effect on the size of the contact patch is also analyzed.
Figure 6 shows the stick-slip distributions of the contact patch for different elastic moduli when the layer thickness is 0.2 mm.The green nodes are in the slip area, and the blue nodes are in the stick area.It can also be observed that the size of the contact patch tends to increase as the elastic modulus increases, although very slightly.Moreover, the shapes of the slip and stick regions are very similar, but the ratio increases from 56.6% to 59.1% when the elastic modulus changes from E 0 to 4E 0 .This means that the ratio of the slip area tends to increase as the layer becomes harder.
To compare the stick-slip conditions and the contact patch sizes for the case with the 0.2 mm thickness, Figure 7 gives the stick-slip distribution of the contact patch for the four different elastic moduli when the layer thickness is 6 mm.The ratio of the slip area gradually reduces from 56.6% for E 0 to 50.0% for 4E 0 .At the same time, the contact patch becomes smaller.This could be because the harder layer produces smaller deformation when the layer becomes thicker, resulting in a smaller slip area and a smaller contact patch.In this case, the contact pressure and surface shear stress become larger, as shown in Figures 4 and 5.In summary, the layer can indeed influence the size of the contact patch and the ratio of the slip and stick.In general, the size of the contact patch is smaller for a harder layer.However, the effect on the ratios of the slip and stick areas depends on the thickness of the layer.

Layer Effects on Subsurface Stresses
This section analyses the layer effects on subsurface stresses.As two key stresses, the von Mises (v-m) stress and the maximum principal stress ( ) are investigated.In particular, the most critical cases-the 0.2 mm and 6 mm thicknesses-are analyzed.Similar to Sections 4.1 and 4.2, the analysis is focused on the stress status when the wheel center is 0.3 m on the x-axis.

Layer Effects on v-m Stress
Figure 8 gives the calculated v-m stresses in the longitudinal cross-section along the central line of the contact patch for different elastic moduli when the thickness is 0.2 mm.If no layer is present, the maximum v-m stress is 665.8MPa, which appears approximately 2 mm below the surface, as shown in Figure 8a.As the elastic modulus increases, the maximum v-m stress also increases.When the elastic modulus is increased to 4 , the v-m stress is maximum, at a value of 1015 MPa-an increase of 52%.In this case, the maximum v-m stress moves to the rail surface.Therefore, an increase in the elastic modulus can produce a larger v-m stress, and it simultaneously tends to move the maximum v-m stress to the surface.This could be because the stress tends to concentrate on the top layer when the elastic modulus becomes larger.In summary, the layer can indeed influence the size of the contact patch and the ratio of the slip and stick.In general, the size of the contact patch is smaller for a harder layer.However, the effect on the ratios the slip and stick areas depends on the thickness of the layer.

Layer Effects on Subsurface Stresses
This section analyses the layer effects on subsurface stresses.As two key stresses, the von Mises (v-m) stress and the maximum principal stress ( ) are investigated.In particular, the most critical cases-the 0.2 mm and 6 mm thicknesses-are analyzed.Similar to Sections 4.1 and 4.2, the analysis is focused on the stress status when the wheel center is 0.3 m on the x-axis.

Layer Effects on v-m Stress
Figure 8 gives the calculated v-m stresses in the longitudinal cross-section along the central line of the contact patch for different elastic moduli when the thickness is 0.2 mm.If no layer is present, the maximum v-m stress is 665.8MPa, which appears approximately 2 mm below the surface, as shown in Figure 8a.As the elastic modulus increases, the maximum v-m stress also increases.When the elastic modulus is increased to 4 , the v-m stress is maximum, at a value of 1015 MPa-an increase of 52%.In this case, the maximum v-m stress moves to the rail surface.Therefore, an increase in the elastic modulus can produce a larger v-m stress, and it simultaneously tends to move the maximum v-m stress to the surface.This could be because the stress tends to concentrate on the top layer when the elastic modulus becomes larger.In summary, the layer can indeed influence the size of the contact patch and the ratio of the slip and stick.In general, the size of the contact patch is smaller for a harder layer.However, the effect on the ratios of the slip and stick areas depends on the thickness of the layer.

Layer Effects on Subsurface Stresses
This section analyses the layer effects on subsurface stresses.As two key stresses, the von Mises (v-m) stress and the maximum principal stress (σ 1 ) are investigated.In particular, the most critical cases-the 0.2 mm and 6 mm thicknesses-are analyzed.Similar to Sections 4.1 and 4.2, the analysis is focused on the stress status when the wheel center is 0.3 m on the x-axis.

Layer Effects on v-m Stress
Figure 8 gives the calculated v-m stresses in the longitudinal cross-section along the central line of the contact patch for different elastic moduli when the thickness is 0.2 mm.If no layer is present, the maximum v-m stress is 665.8MPa, which appears approximately 2 mm below the surface, as shown in Figure 8a.As the elastic modulus increases, the maximum v-m stress also increases.When the elastic modulus is increased to 4E 0 , the v-m stress is maximum, at a value of 1015 MPa-an increase of 52%.In this case, the maximum v-m stress moves to the rail surface.Therefore, an increase in the elastic modulus can produce a larger v-m stress, and it simultaneously tends to move the maximum v-m stress to the surface.This could be because the stress tends to concentrate on the top layer when the elastic modulus becomes larger.To examine if the layer effects on the subsurface stress are associated with the layer thickness, the results of subsurface stresses at 6 mm thickness are also presented, as shown in Figure 9.The maximum v-m stress increases from 665.8 to 751.5 MPa when the elastic modulus changed from  to 4 ; this represents an increase of approximately 13%.This is much smaller than that of the 0.2 mm thickness, illustrated in Figure 8d; however, it is still below the surface, roughly in the middle of the layer.This suggests that the elastic modulus can increase the largest v-m stress slightly but cannot significantly change its location in the case of the 6 mm thickness.Therefore, the layer influences the distribution of the subsurface v-m stress.There is a tendency of a harder layer to induce a larger v-m stress and to change the location of the maximum value, whereas the degree of the layer effects is associated with both the elastic modulus and the thickness of the layer.
To examine if the layer effects on the subsurface stress are associated with the layer thickness, the results of subsurface stresses at 6 mm thickness are also presented, as shown in Figure 9.The maximum v-m stress increases from 665.8 to 751.5 MPa when the elastic modulus changed from E 0 to 4E 0 ; this represents an increase of approximately 13%.This is much smaller than that of the 0.2 mm thickness, illustrated in Figure 8d; however, it is still below the surface, roughly in the middle of the layer.This suggests that the elastic modulus can increase the largest v-m stress slightly but cannot significantly change its location in the case of the 6 mm thickness.To examine if the layer effects on the subsurface stress are associated with the layer thickness, the results of subsurface stresses at 6 mm thickness are also presented, as shown in Figure 9.The maximum v-m stress increases from 665.8 to 751.5 MPa when the elastic modulus changed from  to 4 ; this represents an increase of approximately 13%.This is much smaller than that of the 0.2 mm thickness, illustrated in Figure 8d; however, it is still below the surface, roughly in the middle of the layer.This suggests that the elastic modulus can increase the largest v-m stress slightly but cannot significantly change its location in the case of the 6 mm thickness.Therefore, the layer influences the distribution of the subsurface v-m stress.There is a tendency of a harder layer to induce a larger v-m stress and to change the location of the maximum value, whereas the degree of the layer effects is associated with both the elastic modulus and the thickness of the layer.
Therefore, the layer influences the distribution of the subsurface v-m stress.There is a tendency of a harder layer to induce a larger v-m stress and to change the location of the maximum value, whereas the degree of the layer effects is associated with both the elastic modulus and the thickness of the layer.

Layer Effects on σ 1
Figures 10 and 11 give the calculated σ 1 for the four different elastic moduli when the thicknesses are 0.2 mm and 6 mm, respectively.It is worth noting that the material is in tension if σ 1 is positive; otherwise, the material is in compression.The contact patch is located on the surface during compression.In the case without the layer, the tension occurs outside of the contact patch (represented by the red region in Figure 10a).This phenomenon is in line with that observed in [32].The largest value of tensile σ 1 is 193.8MPa.

Layer Effects on σ
Figures 10 and 11 give the calculated  for the four different elastic moduli when the thicknesses are 0.2 mm and 6 mm, respectively.It is worth noting that the material is in tension if  is positive; otherwise, the material is in compression.The contact patch is located on the surface during compression.In the case without the layer, the tension occurs outside of the contact patch (represented by the red region in Figure 10a).This phenomenon is in line with that observed in [32].The largest value of tensile  is 193.8MPa.
In the case with the 0.2 mm thickness, as the elastic modulus increases, the tensile  becomes larger, as shown in Figure 10.When it increases from  to 4 , the largest tensile  increase from 193.8 MPa to 627.8 MPa is observed, representing an increase of more than 200%.The tension region also tends to move to the topmost surface.Therefore, at a 0.2 mm thickness, a harder layer can induce larger tensile stress and move the tensile region to the surface.The calculated  for the four different elastic moduli when the thickness is 6 mm are presented in Figure 11.In this case, the largest tensile  fluctuates depending on the elastic modulus.The values are 193.8,94.1, 132.9, and 298 MPa for the four elastic moduli, respectively.This means that the influence of the layer may depend on the interaction between the thickness and the elastic modulus.The location of the tensile stress also changes with the elastic modulus.In the case without the layer, the tensile stress is located on the surface (Figure 11a).As the elastic modulus increases to 1.67 and 2 , part of the tensile stress is still located on the surface, while another part appears in the subsurface, specifically below the compression region.When the elastic modulus is increased to 4 , the tensile stresses move completely to the area below the compression region close to the border of the two layers.This could be attributed to the influence of the border.Therefore, the tensile stresses tend to move to the subsurface as the elastic moduli increase; the specific location of the tensile region depends on the elastic modulus.In summary, based on the results of the two layer thicknesses above, the largest tensile stress increases with the elastic modulus for the layer with the 0.2 mm thickness, and the tensile stress is always on the top surface.While the largest tensile stress first decreases and then increases with the elastic modulus for the layer with the 6 mm thickness, the tensile region tends to move to the subsurface area below the compression region.Therefore, the top layer may increase or decrease the magnitude of the tensile stress as the elastic modulus increases, and it may change the location of the tensile area depending on both the elastic modulus and the thickness of the top layer.In the case with the 0.2 mm thickness, as the elastic modulus increases, the tensile σ 1 becomes larger, as shown in Figure 10.When it increases from E 0 to 4E 0 , the largest tensile σ 1 increase from 193.8 MPa to 627.8 MPa is observed, representing an increase of more than 200%.The tension region also tends to move to the topmost surface.Therefore, at a 0.2 mm thickness, a harder layer can induce larger tensile stress and move the tensile region to the surface.
The calculated σ 1 for the four different elastic moduli when the thickness is 6 mm are presented in Figure 11.In this case, the largest tensile σ 1 fluctuates depending on the elastic modulus.The values are 193.8,94.1, 132.9, and 298 MPa for the four elastic moduli, respectively.This means that the influence of the layer may depend on the interaction between the thickness and the elastic modulus.The location of the tensile stress also changes with the elastic modulus.In the case without the layer, the tensile stress is located on the surface (Figure 11a).As the elastic modulus increases to 1.67E 0 and 2E 0 , part of the tensile stress is still located on the surface, while another part appears in the subsurface, specifically below the compression region.When the elastic modulus is increased to 4E 0 , the tensile stresses move completely to the area below the compression region close to the border of the two layers.This could be attributed to the influence of the border.Therefore, the tensile stresses tend to move to the subsurface as the elastic moduli increase; the specific location of the tensile region depends on the elastic modulus.
In summary, based on the results of the two layer thicknesses above, the largest tensile stress increases with the elastic modulus for the layer with the 0.2 mm thickness, and the tensile stress is always on the top surface.While the largest tensile stress first decreases and then increases with the elastic modulus for the layer with the 6 mm thickness, the tensile region tends to move to the subsurface area below the compression region.Therefore, the top layer may increase or decrease the magnitude of the tensile stress as the elastic modulus increases, and it may change the location of the tensile area depending on both the elastic modulus and the thickness of the top layer.

Discussion
Layered structures exist in rails and wheels.The analysis of contact characteristics in previous studies often neglected the layered structures.To determine whether it is necessary to consider the layered effects in actual engineering applications, this study proposed an explicit FE method to study the layer effects.The proposed explicit FE method has high accuracy for wheel-rail rolling contact when an appropriate mesh is applied.The accuracy of frictional rolling has been validated by Zhao and Li [22] for the 0.33, 0.63, and 1.3 mm mesh sizes in the contact surface when the material is homogeneous.Ref. [22] also suggested that the accuracy may be acceptable when the size is 1/10 of the minor axis of the contact patch.To balance the computational cost and accuracy in the present study, the 1 mm mesh is used on the surfaces of the wheel and the rail.This mesh can also give reasonable accuracy compared with Kalker's BEM, as shown in Figure 4.The element size is about 1/10 of the minor axis of the contact patch, which agrees with the conclusion in [22].When the top layer structure is induced, the stresses at the interface change.Kalker's BEM is not appropriate for analyzing cases with layered structures as the accuracy of the contact stresses is difficult to validate directly; however, based on the results in Figures 4 and 5, the harder layer tends to produce larger stresses.This tendency is reasonable, because a harder layer produces a smaller contact patch and, therefore, larger stress.
It is worth noting that there is only one element across the thickness of the layer in the model.To evaluate if this mesh is reasonable for the analysis of the layer effects, two elements across the layer thickness is applied to the case of 1 mm thickness for comparison.The calculated stresses at the interface and the stresses in the subsurface are compared, as shown in Figure 12.As can be seen, the difference is within 1%.Thus, the results of the two mesh types are in good agreement.Moreover, the stresses in the subsurface based on the two mesh types are compared in Figure 13.As can be seen, the stress distribution and the magnitude are very close, and the difference in the maximum value is less than 1%.The small differences indicate that the mesh has a negligible influence on both the stresses at the interface and in the subsurface.Therefore, a mesh with only one element across the thickness is reasonable for studying the layer effects.
based on the two mesh types are compared in Figure 13.As can be seen, the stress distribution and the magnitude are very close, and the difference in the maximum value is less than 1%.The small differences indicate that the mesh has a negligible influence on both the stresses at the interface and in the subsurface.Therefore, a mesh with only one element across the thickness is reasonable for studying the layer effects.When the layer structure appears, the accuracy of the stresses in the subsurface is difficult to validate directly.The present study can be qualitatively validated by comparing studies on multi-layer materials under scratch contact [12,33].In these studies, the contact in scratch is full sliding contact.Some contact characteristics are similar to those in the rolling contact; therefore, they can be used to verify the results in the present study.The authors of [33] studied the stress in multi-layer polymetric systems under scratch contact.They found that a soft substrate was incorporated to induce larger stress fields at the interface of multi-layer materials under scratch contact.The finding is in line with the finding of the effects of elastic modulus of the top layer on the stress at the interface, as shown in Figures 8 and 9. Ref. [12] studied the delamination of multi-layer materials under scratch.They found that tensile maximum principal stress developed at the interface of the multi-layer polymeric structure under the scratch, and the magnitude and direction of the peak tensile maximum principal stress were affected by both the thickness and the material parameters of each layer.Therefore, the findings on the effects of the layer on the maximum principal stress (See Figures 10 and 11) in the present study are in agreement with the observed phenomena in these references [12,33].
The effects of the layer sometimes fluctuate with the elastic modulus, such as the influence on the surface shear stress in Figure 5 and the influence on  in Figure 11.There are two factors that can induce this fluctuation.One is that the effects are determined by both the value of the elastic modulus and the thickness.Another one could be dynamic bution and the magnitude are very close, and the difference in the maximum value is less than 1%.The small differences indicate that the mesh has a negligible influence on both the stresses at the interface and in the subsurface.Therefore, a mesh with only one element across the thickness is reasonable for studying the layer effects.significant as the thickness increases, the trends in the effects become more apparent (see Figure 5e).Therefore, a general trend still exists in the effects.The layer effects should be primarily determined by both the values of the elastic modulus and the thickness.
Based on the numerical results, the influence of the elastic modulus of the top layer is associated with its thickness, such as the influence on the contact pressure in Figure 4 and the influence on σ 1 in Figures 10 and 11.The critical factor behind this could be the ratio between the size of the contact patch and the thickness of the top layer.Therefore, if the size of the contact patch and the thickness are changed, the influence could change.Moreover, the traction coefficient significantly influences the surface shear stress and subsurface stresses.Thus, the layer effects could be also associated with the traction coefficient.In future studies, more traction coefficient values can be considered for comparison.
Although layer effects in rolling contact problems have been extensively studied, e.g., in references [13][14][15][16][17].In these studies, the substrates connected with the layers were assumed to be half-spaces.In reality, the geometries of the wheel and rail are complicated, and the realistic geometries influence the contact characteristics.A robust method that considers realistic geometry is vital in the study of layer effects in engineering applications.Therefore, the explicit FE method was proposed in this study.The realistic geometries of the contact and the substrate were considered.More realistic results were obtained.The results (Figures 4-11) proved that the layer has significantly influence on the contact characteristics, including stresses at the interface, stick-slip behavior and stresses in the subsurface.This implies that the layer structure may need to be considered in more precise practical engineering applications, such as the analysis of wear, RCF and adhesion.For instance, as shown in Figure 10, the layer can induce high tensile stress, which could be the driving force for RCF crack initiation in some conditions.RCF due to WEL could belong to this case, as reported in [34], where crack initiation due to WEL is caused by a brittle fracture.This is why RCF cracks often occur in WEL on rails [35].The present study analyzed a general case of the layer.The precise RCF mechanism still needs detailed analysis to consider more realistic material parameters.Moreover, to study more specific conditions, future studies should include the layer effects of the wheel.
The implication of the results can be addressed from another point of view.At present, some traditional methods are still employed for some wheel-rail contact problems, such as wear prediction (e.g., [36]) and analysis of RCF (e.g., [37]).These methods cannot deal with the problems of layered structures.Since layered structures exist and have a great influence on stress distribution, these methods are not applicable to some practical railway problems.Instead, an FE method should be considered.Precise analyses with consideration of more realistic parameters are still needed to study specific engineering problems.
In the present study, two elastic layers were considered, and the thickness of the layer was assumed.In practice, the material may have multiple layers, and the thickness may vary due to different causes of the layer.Even the material properties gradually change across the thickness.These cases still need to be studied further.In future studies, elastic/plastic material properties should also be considered when dealing with engineering problems.

Conclusions
This study employed an explicit FE method to analyze the elastic layer effects on contact characteristics at the interface and in the subsurface of wheel-rail rolling contact.The top layer structure in rails is considered to have a different elastic modulus from the matrix material.Contact characteristics such as stick-slip behaviors, surface contact stresses, and subsurface stresses are investigated.The accuracy of the solutions with applied mesh was validated by comparing them with Kalker's BEM and other existing studies.The following conclusions can be drawn: The top layer of the rail may alter the contact stresses at the interface.A harder layer induces higher stresses at the interface while reducing the size of the contact patch.The ratios of the stick area and the slip area change as the layer becomes harder.These layer effects tend to become more significant when the thickness of the layer is larger.However, the degrees of these effects are determined by both the value of the elastic modulus and the thickness of the top layer.
The layer also affects the stresses in the subsurface of the rails.There is a tendency for a harder layer to induce larger v-m stress in the subsurface.However, the harder layer may increase or reduce the tensile stress, depending on the thickness.Simultaneously, the layer may change the location of the tensile regions from the surface to the subsurface when the elastic modulus is increased.In general, the degrees of these effects depend on both the elastic modulus and the thickness of the layer.
Layer structures exist in practice, and layer effects are important for wheel-rail rolling contact consequences.The analysis of rolling contact consequences such as wear and RCF may require the inclusion of layer effects by considering more realistic parameters.Moreover, from another point of view, traditional popular methods such as Kalker's BEM and the Hertz theory without consideration of layered structures are not sufficiently accurate for analyzing contact consequences when layered structures exist in the wheel and rail.Instead, methods that consider layered structures should be used.

Figure 1 .
Figure 1.Schematic diagram of wheel-rail rolling contact with a layered structure: (a) Front view; (b) side view.

Figure 1 .
Figure 1.Schematic diagram of wheel-rail rolling contact with a layered structure: (a) Front view; (b) side view.

Lubricants 2023 ,
11, x FOR PEER REVIEW 4 of 16applied to minimize the vibration between the wheel and rail.Its stiffness and damping are 0.95 MN/m and 10.4 kN•s/m, respectively.The weight of the wheel is 660 kg.LMA and CN60 are employed as the wheel and rail profile types, respectively.

Figure 2 .
Figure 2. The mesh of the wheel-rail contact model with a layered structure: (a) Front view; (b) side view.

Figure 2 .
Figure 2. The mesh of the wheel-rail contact model with a layered structure: (a) Front view; (b) side view.

Figure 12 .
Figure 12.Comparison of the contact pressure (P) and surface shear stress (τ).Black lines denote the results of the mesh with one element and green lines denote the results of the mesh with two elements across the thickness.

Figure 13 .
Figure 13.Comparison of v-m stress from two mesh types with a 1 mm thickness: (a) mesh with one element across the thickness of the top layer; (b) mesh with two elements across the thickness of the top layer.

Figure 12 .
Figure 12.Comparison of the contact pressure (P) and surface shear stress (τ).Black lines denote the results of the mesh with one element and green lines denote the results of the mesh with two elements across the thickness.

Figure 12 .
Figure 12.Comparison of the contact pressure (P) and surface shear stress (τ).Black lines denote the results of the mesh with one element and green lines denote the results of the mesh with two elements across the thickness.

Figure 13 .
Figure 13.Comparison of v-m stress from two mesh types with a 1 mm thickness: (a) mesh with one element across the thickness of the top layer; (b) mesh with two elements across the thickness of the top layer.