A simplified model for solving wheel-rail non-Hertzian normal contact problem under the influence of yaw angle

Journal Pre-proof

A simplified model for solving wheel-rail non-Hertzian normal contact problem under the influence of yaw angle Yu Sun , Wanming Zhai , Yunguang Ye , Liming Zhu , Yu Guo PII: DOI: Reference:

S0020-7403(19)33902-5 https://doi.org/10.1016/j.ijmecsci.2020.105554 MS 105554

To appear in:

International Journal of Mechanical Sciences

Received date: Revised date: Accepted date:

14 October 2019 9 February 2020 23 February 2020

Please cite this article as: Yu Sun , Wanming Zhai , Yunguang Ye , Liming Zhu , Yu Guo , A simplified model for solving wheel-rail non-Hertzian normal contact problem under the influence of yaw angle, International Journal of Mechanical Sciences (2020), doi: https://doi.org/10.1016/j.ijmecsci.2020.105554

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.

Highlights 

The approximate expression of the wheel-rail non-Hertzian contact problem considering the yaw angle is proposed.



A simplified wheel-rail non-Hertzian model (MKP-YAW model) that can consider the effect of yaw angle is built.



The MKP-YAW model has a relative good computational accuracy.

A simplified model for solving wheel-rail non-Hertzian normal contact problem under the influence of yaw angle Yu Suna,b, Wanming Zhaib, Yunguang Yec, Liming Zhua and Yu Guod* a

College of Transportation Science & Engineering, Nanjing Tech University, Nanjing 210009, Jiangsu, China

b

State Key Laboratory of Traction Power, Southwest Jiaotong University, Chengdu 610031, China

c

Institute of Land and Sea Transport Systems, Technical University of Berlin, Berlin 10587, Germany

d

Transportation Institute, Inner Mongolia University, Hohhot 010023, Inner Mongolia, China

Abstract: Wheelset yaw motion significantly affects wheel-rail non-Hertzian contact behaviour, but currently, few wheel-rail normal contact models consider this issue. In this paper, a simplified wheel-rail contact model is built to study the wheel-rail non-Hertzian normal contact problem considering the yaw angle. First, three basic assumptions for the approximate calculation of this problem are proposed and are validated by using results calculated by Kalker's variational method (KVM), and the approximate expression of the wheel-rail non-Hertzian contact problem considering the yaw angle is obtained based on these assumptions. Then, the approximate expression is applied to the MKP (Modified Kik–Piotrowski) contact model, and an improved MKP model that considers the wheelset yaw angle is proposed, i.e., the MKP-YAW model. Compared with the results calculated by KVM, the computational accuracy of the proposed MKP-YAW model in calculating the wheel-rail contact considering the yaw angle is similar to that of the MKP contact model in calculating the wheel-rail contact without considering the yaw angle. Keywords: Wheel-rail contact, non-Hertzian contact, yaw angle, Kalker's variational method, MKP contact model

1. Introduction Wheel-rail contact simulation is a crucial issue in vehicle system dynamics [1], train-track-bridge interactions [2] and wheel-rail damage analysis [3]. The wheel-rail normal contact model originated from Hertzian contact theory [4,5], i.e., the Hertzian contact model, which is currently the most widely used model in wheel-rail normal contact analysis [6]. However, wheel-rail contact is essentially a non-Hertzian contact problem as the curvatures of the wheel and rail profiles are constantly changing, especially for wheel-rail contact on seriously worn profiles. The non-Hertzian contact model, therefore, may be more suitable for describing the contact behaviour between the wheel and the rail. Due to the complexity of wheel and rail profiles, it is impossible to analytically solve the wheel-rail non-Hertzian contact problem in the same manner as that for Hertzian contact. To address this issue, several numerical calculation models have been proposed, including the finite element method (FEM) [7-11] and boundary element method (BEM) [12-17]. Among them, the FEM can provide the best computational accuracy when solving the wheel-rail non-Hertzian contact problem. The FEM is suitable for solving extremely complicated wheel-rail contact problems, such as conformal contact [17] and elastic-plastic contact [11]. This method, however, is hard to directly apply to vehicle system dynamics simulations due to the very large computational effort required. In terms of the BEM, the Green functions that relate the displacement and the stress of the contact surface are crucial in the calculations [18]. Wheel-rail materials *

Corresponding author. E-mail address: [email protected] (Yu Guo)

are usually assumed to be homogeneous, isotropic, and linear elastic; thus, the elastic half-space assumption is widely adopted in wheel-rail contact analysis [13]. For the elastic half-space, the Green function is given by the Boussinesq-Cerruti formula [5,12,13]. Kalker's variational method (KVM) [13], which is based on the Boussinesq-Cerruti formula, is the most famous BEM in the analysis of wheel-rail contact. KVM is more efficient than the FEM, but this method still cannot meet the requirements of online calculations, even though a solver for accelerating the calculation has been proposed by Vollebregt [14]. The elastic half-space assumption is not satisfied for wheel-rail conformal contact as the contact patch is curved. A description of the wheel-rail conformal contact using BEM can be found in Refs. [15,16]. To overcome the problem of low efficiency in wheel-rail non-Hertzian contact analysis, many approximate methods have been proposed. Specifically, a series of wheel-rail non-Hertzian contact models based on the virtual penetration method have been developed in the past two decades, including the Linder method [19], the Kik–Piotrowski (KP) model [20,21], the Extended Kik–Piotrowski (EKP) model [22], the Modified Kik–Piotrowski (MKP) model [23-25], the Ayasse-Chollet model (STRIPES) [21,26], and the ANALYN model [27]. The computational robustness of the wheel-rail contact model is critical in vehicle system dynamics simulations. The local curvatures of the wheel-rail profiles are adopted in the methods of STRIPES, ANALYN, and EKP, which adversely affects the computational robustness as the local curvatures are usually discontinuous and sometimes fluctuate significantly. The abnormal curvatures need to be smoothed before the wheel-rail contact analysis. However, different smoothing methods often lead to different calculation results, and the smoothing methods are often very subjective, resulting in poor robustness [23]. In the KP and MKP contact models, the local curvatures of wheel-rail profiles are not considered, which results in a higher robustness in calculating the wheel-rail normal force. However, ignoring the influence of the profile curvature on the wheel-rail normal contact could cause worse estimations of the shape of the wheel-rail contact patch and the contact pressure [22,27]. As wheelset yaw motion is inevitable in vehicle system dynamics, this motion will affect the wheel-rail normal contact behaviour. To improve the accuracy of the wheel-rail contact model, it is of critical importance to consider the influence of the yaw angle in the wheel-rail non-Hertzian contact model, especially when the train is running on a small radius curve since a large yaw angle will usually occur. However, currently, few studies have paid attention to this issue. Most of the existing wheel-rail non-Hertzian contact models have difficulty directly calculating the wheel-rail non-Hertzian contact with the yaw angle. Vo et al. [28] studied the wheel-rail contact with different yaw angles by FEM. An et al. [29] proposed a wheel-rail spatial geometric gap estimation method for calculating wheel-rail non-Hertzian contact with the yaw angle by KVM. Liu et al. [22] extended the KP contact model to obtain the EKP contact model, which could consider the influence of the yaw angle. Note that the spatial geometric gap between the wheel-rail contact surface is considered in these models, and the computational effort of the spatial geometric gap is very large for the vehicle system dynamics analysis. Based on the above consideration, it is of great interest to develop an efficient wheel-rail non-Hertzian contact model that is able to consider the wheelset yaw angle. In this paper, a simplified model for solving the wheel-rail normal contact problem that considers the yaw angle is presented. Based on the calculated results of KVM (implemented in a program named ‗S-CONTACT‘ coded by the authors), three basic assumptions relating to the wheel-rail non-Hertzian contact with the yaw angle are proposed and validated. Combining these three basic assumptions together, we obtain the approximate expression for the wheel-rail non-Hertzian contact problem with the yaw angle, which can be further applied to the approximate wheel-rail non-Hertzian model with the yaw angle. The MKP contact model is applied to illuminate the usage of the approximate expression in the approximate wheel-rail non-Hertzian model. The MKP contact

model is combined with the approximate expression to obtain a simplified wheel-rail non-Hertzian model that can consider the effect of the yaw angle, i.e., the MKP-YAW contact model. Two-dimension geometric gap is adopted in the MKP-YAW contact model, which can reduce the computational complexity and effort of the wheel-rail contact geometry calculation. The computational accuracy of the MKP-YAW model is compared with that of the MKP model by taking the results estimated by KVM (S-CONTACT) as the reference.

2. Basic assumptions of wheel-rail non-Hertzian contact with yaw angle To simplify the calculation of the wheel-rail non-Hertzian contact under the influence of the wheelset yaw angle, three basic assumptions are proposed and validated in this section based on the characteristics of wheel-rail contact obtained by KVM, and the approximate expression for the wheel-rail non-Hertzian contact problem with yaw angle is obtained.

2.1 Introductory remarks The coordinate system O-xyz, named the wheel-rail contact coordinates, is defined in Fig. 1. The x-axis points in the rolling direction, the z-axis is normal to the wheel-rail contact plane, and the y-axis is obtained from the right-hand rule, where the coordinate origin O is the initial contact point. The following wheel-rail normal contact analysis is established in the wheel-rail contact coordinates O-xyz.

x y O

z Figure 1. Definition of wheel-rail contact coordinate system. x

(b) (a)

Pressure at the main axis xp(y)

Contact patch

xr(y)

ar(y) pz(x,y) x y

y

xp(y)

O yl

yr al(y)

Main axis of the normal pressure

xl(y)

Figure 2. Wheel-rail normal contact with the yaw angle. (a) Distribution of wheel-rail normal contact pressure, (b) shape of the wheel-rail contact patch.

The contact patch and the normal pressure are always symmetric about the y-axis when the yaw angle is ignored. However, this conclusion cannot be applied to the wheel-rail contact with the yaw angle. A typical case of wheel-rail normal contact with the yaw angle is shown in Fig. 2. The normal contact pressure and the contact patch are distorted due to the influence of the yaw angle. The maximum contact pressures of different lateral cross-sections (cross-sections parallel to the x-axis and with different y-coordinates) are no longer located at the main contour plane x=0. Define the curve inside the contact patch where the corresponding normal pressures are the largest for different lateral cross-sections as the "main axis of the normal distance", expressed as xp(y). The main axis of the normal pressure xp(y) satisfies

pz ( xp y , y)

p ) , z ( x, y

x,y

C

(1)

where C represents the area of the contact patch. As shown in Fig. 2(b), the left and right edges of the contact patch are defined as yl and yr, the rear and front edges of the contact patch are defined as xl(y) and xr(y), and the distances between the rear and front edges of the contact patch and the main axis of the normal pressure are defined as al(y) and ar(y), respectively. Therefore, al(y) and ar(y) can be expressed as

al ,r  y   xl ,r  y   x p  y 

(2)

2.2 Design of the simulation conditions The effects of the yaw angle on the shape of the contact patch and the normal contact pressure pz(x,y) are discussed in detail in Section 2.3, and the simulation conditions are first described in this section. KVM is a boundary element method for solving elastic half-space contact problems. Since KVM has a high computational accuracy for wheel-rail non-Hertzian contact analysis, the self-designed program named S-CONTACT based on KVM [13,23] is adopted to study the basic characteristics of the wheel-rail non-Hertzian contact with the yaw angle. The rectangular element and fixed normal pressure inside each element are adopted in KVM. The element size used in the latter simulation is 0.2 mm in both the longitudinal and lateral directions. The material properties of the wheel and the rail are assumed to be the same, the shear modulus is 82 GPa and Poisson's ratio is 0.28. 30 1st type profile 2nd type profile 3rd type profile

25 2

Z (mm)

20

1 0

15

-1

10

-2 -20

-10

0

10

20

5 0 -5 -60

-40

-20

0 Y (mm)

20

40

60

Figure 3. Comparison of wheel profile used in simulation. 1st type profile: standard wheel profile LMA, 2nd type profile: hollow-worn LMA with a running mileage of 83,000 km, 3rd type profile: hollow-worn LMA with a running mileage of 192,000 km.

To reflect the influence of the complex wheel-rail contact states, the contact behaviour between the standard Chinese rail profile CHN 60 and three types of wheel profiles, including one standard wheel profile LMA and two worn LMA profiles, are investigated. LMA is one of the most widely used wheel profile in

China, and the two worn wheel profiles are measured from high-speed EMUs with running mileages of 83,000 km and 192,000 km [30]. The standard LMA profile and the worn LMA profiles with running mileages of 83,000 km and 192,000 km are named the first type, second type, and third type of wheel profiles, respectively. A comparison of these three profiles is shown in Fig. 3. The profiles show that the wear region is primarily located at the wheel tread around the nominal rolling circle. The maximum wear depth is 0.33 mm for the second type profile and 0.72 mm for the third type profile. These three profiles were also applied to study the wheel-rail normal contact without yaw angle in Ref. [23]. Wheel-rail normal contact with a fixed approach under the conditions of different wheelset lateral displacements and yaw angles are adopted in the analysis. The specific wheel-rail approach of 0.08 mm is selected. The wheelset lateral displacement ranges from -6 mm to 6 mm at an interval of 0.05 mm. According to Ref. [22], wheel-rail contact with a yaw angle of 0.05 rad can reflect the operating conditions when the vehicle runs in sharp curves. For this reason, the selected yaw angle ranges from 0.005 rad to 0.05 rad at an interval of 0.005 rad. In conclusion, for each of the three wheel-rail contact conditions, there are 241 wheelset lateral displacement cases and 10 wheelset yaw angle cases, with 2410 combined calculation cases in total. According to the principle of symmetry, the wheel-rail contact with a negative yaw angle can be expressed by the wheel-rail contact with a positive yaw angle. The equivalence relationships are as follows: (1) The wheel-rail normal pressure satisfies

pz ( x, y, Yw ,  w )  pz ( x, y, Yw , w )

(3)

(2) The shape of the wheel-rail contact patch satisfies

xp  y,Yw ,

  xp  y Y, w, w

(4)

w

al  y, Yw ,  w   ar  y, Yw , w 

(5)

where Yw represents the wheelset lateral displacement and w represents the wheelset yaw angle. Considering the equivalence relationship given above, wheel-rail contact with a negative yaw angle can be ignored in the analysis.

2.3 Basic assumptions of wheel-rail non-Hertzian contact under the effect of yaw angle The effect of the yaw angle on the wheel-rail normal contact is limited since the yaw angle is relatively small. For this reason, a series of assumptions can be applied in the simplified calculation of the wheel-rail normal contact with the yaw angle according to the conclusions of normal contact without the yaw angle. Three assumptions for the wheel-rail non-Hertzian contact solution are proposed below, and the results obtained by KVM under the simulation conditions described in Section 2.2 are applied to verify the assumptions. 2.3.1 Assumption 2.1 Define the curve inside the contact patch at which the normal distance between undeformed wheel-rail surfaces is the shortest for different lateral cross-sections as the "main axis of the normal distance", expressed as xd(y). The main axis of the normal distance xd(y) satisfies zw r( xd  y , y)

z ( x, y ) ,

w r

x,y

where zwr is the normal distance between the undeformed wheel-rail surfaces.

C

(6)

Assumption 2.1: The main axis of the normal pressure coincides with the main axis of the normal distance. Assumption 2.1 can also be expressed as

xp ( y) xd ( y) , ly

 y

r

(7)

y

The assumption that the rail is straight and smooth in the rolling direction is adopted; this assumption is widely used to calculate the wheel-rail contact geometry. Based on this assumption, the normal vector of the wheel surface at the initial wheel-rail contact point should be parallel to axis O0z0 in Fig. 1. All the points satisfy the requirements of the normal vector belonging to a wheel-rail potential contact space curve on the wheel surface. This space curve is named the “contact locus” [31]. Note that the points satisfying the requirements of the normal vector are equivalent to those satisfying the rule that the normal distance between the undeformed wheel-rail surfaces is the shortest of the corresponding lateral cross-section. Therefore, the main axis of the normal distance xd(y) can be obtained by solving the expression of the contact locus. The contact locus can be analytically solved by Wang‘s method [31-33]. Figure 4 sketches the wheel-rail contact geometry calculation based on Wang‘s method. Four coordinate systems can be seen in Fig. 4, which include the absolute coordinates (O0-x0y0z0), the transitional wheel coordinates (O1-x1‟y1‟z1‟), the wheel coordinates (O1-x1y1z1), and the wheel-rail contact coordinates (O-xyz). For the absolute coordinates, the coordinate origin O0 is located at the track centre, the x0-axis points in the rolling direction, the z0-axis points to vertically downward, and the lateral y0-axis points along the direction of the track transverse obtained from the right-hand rule. The transitional wheel coordinates (O1-x1‟y1‟z1‟) are parallel to the coordinates O0-x0y0z0, and the coordinate origin O1 is located at the centre of the wheelset. The coordinates of the wheelset centre O1 in the absolute coordinates are (0, Yw, Zw), where Zw is the vertical displacement of the wheelset. The wheel coordinates (O1-x1y1z1) are obtained by rotating w (wheelset yaw angle) along the O1-z1‟-axis and rotating w (wheelset rotation angle) in the O1-x1‟-axis by the coordinates O1-x1‟y1‟z1‟. The lateral y1-axis is parallel to the axle of the wheelset. The wheel-rail contact coordinates (O-xyz) are introduced in Section 2.1. In wheel-rail contact analysis, the wheel profile is always described by the coordinates (yw, Rw, w), where yw is the lateral coordinate of the wheel profile in coordinates O1-x1y1z1, Rw is the corresponding local rolling radius, and w is the corresponding local contact angle of the wheel surface. Wheel rolling circle

x0 x1'

w x1

Rw

Wheelset central line

x

B y0

O0

Zw Yw

y1

O

y

w w

O1

Rail surface

y1'

z z0

z1'

z1

Figure 4. Schematic diagram of wheel-rail contact geometry solved by Wang‘s method.

According to Wang‘s method, in the absolute coordinates (O0-x0y0z0), the coordinate of the contact locus (xd0, yd0, zd0) can be expressed as

 xd 0  y w  y lwx  lx R wy  wtan   yw  w   y  y   y l  Y  Rw  yw  l 2l tan   y   l 1  l 2 (1  tan 2   y )   d0 w w y w x y w w z x w w  1  lx2     z  y   y l  Z  Rw  yw  l 2l tan   y   l 1  l 2 (1  tan 2   y )  d0 w w z w x z w w y x w w   1  lx2  

(8)

where lx, ly, and lz are the direction vectors of the wheel coordinates related to the absolute coordinates, which can be written as  lx   cos w sin w   l y  cos w cos w l  sin  w z

(9)

The wheel-rail initial contact point can be obtained by finding the minimum vertical distance between the contact locus and the rail profile in the absolute coordinates. Descriptions of the wheel-rail contact geometry can be found in Ref. [32,34] in detail. The lateral coordinate of the initial contact point is defined as „yc‟ in the wheel coordinates. According to the conversion relationship between the absolute coordinates (O0-x0y0z0) and the wheel-rail contact coordinates (O-xyz), the main axis of the normal distance can be expressed as  xd  yw   xd 0  yw   xd 0  yc    yd  yw    yd 0  yw   yd 0  yc  / cos  w  w 

(10)

The main axis of the normal pressure can be obtained by combining Assumption 2.1 (Eq. (7)) and the analytical solution of the main axis of the normal distance (Eqs. (8)-(10)). Element central

xd(y) xp(y)

x (mm)

x (mm)

xd(y)±xe

xd(y)+xe xd(y) xd(y)-xe xp(y)

Contact area

y (mm)

y (mm)

Figure 5. Comparison between xp(y) and xd(y) under the simulation condition of the first type of wheel, Yw=5 mm, and

w=0.05 rad.

Figure 5 shows a comparison of the results for the main axis of the normal pressure xp(y) and the main axis of the normal distance xd(y), in which the wheel-rail normal contact between the first type of wheel (LMA) and the rail CHN 60 with a lateral displacement Yw=5 mm and yaw angle w=0.05 rad is adopted in the analysis. xp(y) is obtained by KVM, and xd(y) is obtained by Wang‘s method. Since the rectangular element with a fixed size is used in KVM, xp(y) appears as a jagged line. The discrete mesh used in KVM causes a certain computational error in the calculation of xp(y), and the calculation error induced by the discrete mesh is basically within one element size in the rolling direction ( xe). As shown in Fig. 5, xd(y) is relatively in good accordance with xp(y). xp(y) is distributed on both sides of xd(y), which is the central axis, and the calculation error between xp(y) and xd(y) is basically within xe. Considering the calculation error induced by the discrete mesh, the normalized calculation error between xp(y) and xd(y) is shown here to validate Assumption 2.1, which can be expressed by

ERA1 ( y) 

xd ( y)  x p ( y)

 xe

, yl  y  yr

(11)

Wheel-rail contact between the three types of wheels and rail CHN 60 with different values of wheelset lateral displacement and yaw angle are calculated by KVM and Wang‘s method. The detailed description concerning the wheel-rail contact conditions can be seen in Section 2.2. For each type of wheel, the normalized calculation errors of all the KVM lateral meshes inside the contact patch for the Yw & w combined cases (2410 cases) are counted. Figure 6 shows the probability distribution of the normalized calculation error of Assumption 2.1. As shown in Fig. 6(a), the calculation errors between xp(y) and xd(y) are very small for all three types of wheel cases. The calculation errors are mainly within xe, and almost all the calculation errors are less than 3∗xe. Moreover, the probability distribution of the normalized calculation error is almost symmetrical about ‗ERA1=0‘, which indicates that the calculation errors may be caused by the discrete mesh of KVM. The calculation error increases with an increasing degree of wheel wear. The probabilities of a calculation error within xe are 95.7%, 84.1%, and 79.2% for the first, second, and third types of wheel profiles, respectively. (b)

(a) 1.2

Probability density

1.0

Cumulative probability density

1st type profile 2nd type profile 3rd type profile

0.8 0.6 0.4 0.2 0.0 -3

-2

-1

0 ERA1

1

2

3

1.0 0.8 0.6 0.4 1st type profile 2nd type profile 3rd type profile

0.2 0.0 0.0

0.5

1.0

1.5 |ERA1|

2.0

2.5

3.0

Figure 6. Normalized calculation error of Assumption 2.1. (a) Probability density function of ERA1, (b) cumulative probability density function of the absolute value of ERA1.

As stated above, Assumption 2.1 can be adopted in the simplified calculation of wheel-rail non-Hertzian contact under the influence of the yaw angle. 2.3.2 Assumption 2.2 Remark 2.1: Based on Assumption 2.1, al(y) and ar(y) denote the distances between the rear and front edges of the contact patch and the main axis of the normal distance xd(y), respectively, in the following analysis. Eq. (2) can be written as al ,r  y   xl ,r  y   xd  y 

(12)

The equation “al(y)=ar(y)” is obviously correct for the wheel-rail contact without the yaw angle. Assuming that the wheelset yaw motion has little effect on the relationship between al(y) and ar(y) since the yaw angle is relatively small, we can make the following assumption: Assumption 2.2: The distances between the rear and front edges of the contact patch and the main axis of the normal distance are the same, namely, al(y)=ar(y). Figure 7 shows a comparison of al(y) and ar(y) obtained by KVM under simulation condition that is the same as that adopted in Fig. 5 (first type of wheel contact with CHN 60, Yw=5 mm, andw=0.05). As shown in Fig. 7, al(y) and ar(y) are in good agreement at different lateral positions. Assumption 2.2 could be applied

in this contact situation with high computational accuracy. 6 al

al,r (mm)

5

ar

4 3 2 1 0 -8

-4

0 y (mm)

4

8

Figure 7. Comparison of al(y) and ar(y) under the simulation condition of the first type of wheel, Yw=5 mm, and w=0.05 rad.

To comprehensively validate Assumption 2.2, the normalized calculation error between al(y) and ar(y) is introduced, which can be expressed by a ( y )  ar ( y ) ERA2 ( y )  l , yl  y  yr (13)  xe Under the simulation conditions adopted in the calculation of Fig. 6, the probability distribution of the normalized calculation errors of Assumption 2.2 is shown in Fig. 8. Even though the calculation errors of Assumption 2.2 are relatively larger than those of Assumption 2.1 shown in Fig. 6, the errors are still very small and acceptable in the simplified calculation. The probabilities of calculation error between al(y) and ar(y) within xe are 74.3%, 67.5%, and 61.7% for the first, second, and third types of wheels, respectively. The calculation errors between al(y) and ar(y) less than 3∗xe satisfy most of the simulation conditions. Fig. 8(a) also shows the approximate symmetric distribution of the probability distribution of ERA2 with the symmetry axis at ‗ERA2=0‘. (b)

Cumulative probability density

(a) 0.5

Probability density

0.4 0.3 0.2 0.1 1st type profile 2nd type profile 3rd type profile

0.0 -0.1 -3

-2

-1

0

1

2

3

1.0 0.8 0.6 0.4

1st type profile 2nd type profile 3rd type profile

0.2 0.0 0.0

ERA2 (%)

0.5

1.0

1.5 2.0 |ERA2| (%)

2.5

3.0

Figure 8. Normalized calculation error of Assumption 2.1. (a) Probability density function of ERA2 and (b) cumulative probability density function of the absolute value of ERA2.

Remark 2.2: Assumption 2.2 is adopted in the following analysis, and we define a  y   al  y   ar  y 

(14)

2.3.3 Assumption 2.3 Since the yaw angle is relatively small, the wheelset yaw motion has little effect on the expression of the normal pressure distribution along the rolling direction, and the following assumption can be obtained: Assumption 2.3: The normal pressure distribution is semi-elliptical in the rolling direction.

The accuracy of Assumption 2.3 is evaluated by combining the assumptions proposed previously. Taking Assumption 2.1, Assumption 2.2, and Assumption 2.3 all into consideration, the approximate expression of the wheel-rail normal contact problem under the influence of the yaw angle can be reduced as: (1) wheel-rail normal pressure

pz  x, y   pz  xd ( y), y 

 x  xd ( y)  1    a( y ) 

2

(15a)

(2) wheel-rail contact area C a  y  , yl  y  yr

C : xl , r ( y )  xd ( y )

(15b)

To evaluate the accuracy of the approximate expression (Eq. (15)), the wheel-rail non-Hertzian contact calculated by KVM is equivalent to the following expressions based on Eq. (15) pzSK  x, y   pzK  xd ( y ), y 

 x  xd ( y )  1    aSK ( y ) 

2

(16a)

aSK ( y )   alK ( y )  arK ( y )  / 2

(16b)

where the subscript „K‘ denotes the variable calculated by KVM and the subscript „SK‘ denotes the variable obtained by the equivalent expression (Eq. (16)). The contact patch and the contact pressure obtained by the equivalent expressions are compared with those calculated by KVM. The function of the ―overlap ratio‖ is introduced to evaluate the computational accuracy. The definition of the function of the overlap ratio is described as follows. Remark 2.3: The overlap ratio of two closed spaces C1 and C2 is denoted as twice the ratio of the area (two-dimensional space) or volume (three-dimensional space) of space C1∩C2 to the total area/volume of C1 and C2, namely, ORC1 / C2 

2  D  C1

C2 

(17)

D  C1   D  C2 

where D(C) represents the area of the closed two-dimensional space C or the volume of the closed three-dimensional space C. The overlap ratio can be used to reflect the agreement of the position and shape of two closed spaces. According to the definition of the overlap ratio, the overlap ratio of the wheel-rail contact patch obtained by KVM and the equivalent expression can be derived by ORCSK / C K 



2  D  CSK

CK 

D  CSK   D  CK  2 

yrK yK

min  x ( y)  a d

  min  x yrK



ylK

d

yrK ylK

SK

 y  , xrK ( y )   max  xd ( y)  aSK  y  , xlK ( y )  dy

2  aSK  y  dy  

yrK ylK

 alK ( y )  arK ( y )  dy

(18)



( y )  aSK  y  , xrK ( y )   max  xd ( y )  aSK  y  , xlK ( y )  dy 2 

yrK ylK

aSK  y  dy

where CK and CSK represent the contact area obtained by KVM and the equivalent expression, respectively. The overlap ratio of the wheel-rail contact pressure obtained by KVM and the equivalent expression can

be derived by

ORVSK /V K 

2  D VSK

VSK 

D VSK   D VK 



2 

yrK ylK



min  xd ( y )  aSK  y , xrK ( y ) 

max  xd ( y )  aSK  y , xlK ( y ) 

yrK

xd ( y )  aSK  y 

ylK

xd ( y )  aSK  y 

 

min  pzSK  x, y  , pzK  x, y   dxdy

pzSK  x, y  dxdy  

yrK ylK



xrK ( y )

xlK ( y )

pzK  x, y  dxdy

(19)

where VK and VSK represent the three-dimensional space formed by the plane z=0 and the contact pressure distribution calculated by KVM and the equivalent expression, respectively. Figure 9 shows the overlap ratios of the wheel-rail contact patch and the contact pressure calculated by KVM and the equivalent expression (Eq. (16)) for different wheel types, wheelset lateral displacements and yaw angle cases. The overlap ratios are extremely high for all the simulation cases. Overall, the overlap ratio decreases slightly with increasing yaw angle, and the overlap ratio of the standard wheel (first type of wheel) is slightly better than those of the worn wheels (second and third types of wheels). The lowest overlap ratio of the wheel-rail contact patch is 96.45%, which occurs for the second type of wheel, Yw=3.55 mm, and w=0.05 rad. The lowest overlap ratio for the wheel-rail contact pressure is 97.75%, which also occurs for the second type of wheel, Yw=5.05 mm, and w=0.05 rad. Figure 10 shows the comparison of the contact area for the lowest overlap ratio case of the wheel-rail contact patch. As shown in Fig. 10, the contact area calculated by the equivalent expression (CSK) still has good agreement with that calculated by KVM (CK), even though the overlap ratio is the lowest in all contact conditions.

(a1)

(a2)

(b1)

(b2)

(c1)

(c2)

Figure 9. Overlap ratios of wheel-rail contact patch (left) and contact pressure (right) obtained by KVM and the equivalent expression. (a) First type of wheel, (b) second type of wheel, (c) third type of wheel. 8 CK CSK

x (mm)

4

CK∩CSK

0

-4 ORCSK/CK= 96.45%

-8 -4

-2

0

2 y (mm)

4

6

8

Figure 10. Comparison of the contact areas calculated by KVM (CK) and the equivalent expression (CSK) under the simulation condition of the second type of wheel, Yw=3.55 mm, and w=0.05 rad.

Consequently, the approximate expression for the wheel-rail non-Hertzian contact problem (Eq. (15)) considering the yaw angle can be adopted in the simplified calculation with a high accuracy. Based on the approximate expression, the wheel-rail non-Hertzian contact with yaw angle problem is simplified to solve the shape of the contact patch (yl,r, a(y) in Eq. (15a)) and the contact pressure of the main axis of the normal distance (pz(xd(y),y) in Eq. (15b)).

3. Application of the approximate expression in simplified wheel-rail non-Hertzian contact calculation The approximate expression (Eq. (15)) based on the basic assumptions is applied in the wheel-rail non-Hertzian model. A robust non-Hertzian model proposed in the authors‘ previous work [23-25], which is named the MKP contact model, is combined with the approximate expression to obtain the improved MKP contact model, which can consider the influence of the yaw angle. The improved MKP contact model is abbreviated as the MKP-YAW contact model.

3.1. MKP contact model The MKP contact model is a simplified wheel-rail non-Hertzian contact model that is based on the virtual penetration method. This model improves the computational accuracy of the original KP model and exhibits better computational accuracy relative to the KP model. According to the MKP contact model, the wheel-rail normal contact pressure can be written as

p z ( x, y ) 

E  wr  y  2 1 2





a2  y   x2 yl

a  

yr a 

(20)

a 2     2

 2    y 

2

d d

where E and  are the elastic modulus and Poisson's ratio of the wheel/rail, respectively, and wr(y) is the elastic deformation of the coordinates (0, y) [23].

The wheel-rail normal contact force can be obtained by integrating the normal contact pressure inside the contact patch, namely, N

 2E 4 1 2





a 2  y   wr  y 

yl yr

yl

yr a 

 2    y 

(21)

dy

a 2     2

a  

2

d  d

In the case of multi-point contact, the wheel-rail normal force can be expressed as N

 2E 4 1 2





 i 1 y Nc

y

a 2  y   wr  y 

li

ri

a  

 j 1 y lj a  Nc

y

rj

dy

a 2     2

 2    y 

2

(22)

d  d

where the subscripts „i‟ and ‟j‟ represent the serial number of the contact patch and Nc is the total number of contact patches. A more detailed description concerning the MKP contact model can be found in Refs. [23-25].

3.2. MKP-YAW contact model The approximate expression of the wheel-rail non-Hertzian contact problem (Eq. (15)) is applied to the MKP-YAW contact model, which allows the MKP model to calculate the wheel-rail non-Hertzian contact with the yaw angle. 3.2.1 Estimation of the wheel-rail contact patch in the MKP-YAW contact model In the MKP-YAW contact model, the wheel-rail contact is estimated by almost the same method as that in the MKP model. The main axis of the normal distance, which is regarded as the centre of the contact patch of the approximate expression, is used in the MKP-YAW contact model to replace the axis of symmetry of the contact patch (x=0) of the MKP model. Assuming that the shape of the contact patch can be determined by the virtual interpenetration area [20,23], the left or right edge of the original contact patch 𝑦̅𝑙,𝑟 satisfies





zwr xd  yl ,r  , yl , r  

(23)

where  is the reduction factor, which is selected to be 0.5 in the MKP-YAW contact model. The half longitudinal length of the original contact patch satisfies

a  y   2R  y    zwr  xd  y  , y  , yl  y  yr

(24)

The contact patch is corrected by Hertzian theory, and the corrected contact patch can be expressed as

a ( y )  a y( )y,l  y  yr

(25)

where  and  are the ratios of the corrected and the original size of the contact patch, respectively see the Appendix. To ensure the compressed state of the contact area, the contact patch needs to be further modified, and a detailed description of the correctional method can be found in Ref. [23]. 3.2.2 Extension of the KP model Using the assumptions adopted in the KP contact model, the wheel-rail contact pressure at the main axis of the normal pressure (or the main axis of the normal distance) is assumed to be similarly distributed as the

half longitudinal length of the contact patch, namely, p ( x (0),0) pz ( xd ( y ), y )= z d a y a  0

(26)

Combined with Eq. (15), the normal pressure of the KP model with yaw angle can be expressed as  x  xd  y   p ( x (0),0) pz ( x, y )= z d a y 1   a  0  a  y  

2

(27)

The wheel and rail are usually simplified into the elastic half-space in the wheel-rail contact analysis. Based on the elastic half-space assumption, the Boussinesq-Cerruti formula can be applied to describe the relationship between the surface pressure and the surface elastic deformation. In the KP model, the relationship between the wheel-rail normal pressure distribution (Eq. (27)) and the elastic deformation of the initial contact point is adopted. Consequently, the simulation result is sensitively affected by the initial contact point [23]. As an extension of the KP method, the influence of the wheel-rail normal pressure (Eq. (27)) on the elastic deformation of the main axis of the normal distance is included. When the elastic deformation of the point on the main axis of the normal distance (xd(y‟), y‟) and the wheel-rail normal pressure distribution inside the contact patch satisfy the Boussinesq-Cerruti formula [5,12,13], the corresponding normal pressure of the initial contact point is defined as P0(y‟), and the normal pressure distribution is defined as Pz (x, y, y‟). According to Eq. (27), Pz (x, y, y‟) can be expressed as Pz ( x, y, y )  P0  y  

a y a  0

 x  xd  y   1    a  y  

2

(28)

Based on the Boussinesq-Cerruti formula, the total elastic deformation of the wheel-rail surface at the position of (xd(y‟), y‟) can be expressed as uwr  xd  y   , y     wr  y   









2 1   2 P0  y  



 Ea  0 



2 1   2 P0  y  

 Ea  0 



2 1 2

E yl

x



   a  

yr xd  a  yl

d

a  

yr a 

yl

x

d

P  , , y  

   a  

yr xd  a 

2   xd  y     y2

a 2      xd   

2

2   xd  y     y2

xl2     2

d  d

d  d

2   xd    xd  y     y2

(29)

d  d

where uwr is the elastic deformation of the wheel-rail contact surface and wr (y‟) is the elastic deformation of the point (xd(y‟), y‟) on the main axis of the normal distance, which can be obtained by the wheel-rail contact geometry

 wr  y    zwr  xd  y , y

(30)

where the normal distance between the undeformed wheel-rail surfaces zwr can be written as zwr  xd  y  , y    zw  xd  y  , y   zr  xd  y  , y  cos  w  w    z w  xd  y  , y   zr  y   cos  w  w  (31)

where zw(xd(y),y) describes the position of the wheel profile on the main axis of the normal distance, which can be obtained by the coordinates of the contact locus, and zr(y) describes the position of the rail profile [32].

According to Eq. (29), the expression of P0(y‟) is deduced to be P0  y   =

 Ea  0   wr  y  



2 1 2



/

yl yr

a 2     2

a  

a 

  xd    xd  y       y   2

2

d  d

(32)

3.2.3 Normal contact force of the MKP-YAW contact model Analogous to the assumptions used in the MKP model, the wheel-rail normal contact pressure under the influence of the yaw angle can be expressed as pz ( x, y )  Pz ( x, y, y ) 

a 2  y    x  xd  y  

 E wr  y 



2 1 2



yl

(33)

a 2     2

a  

yr a 

2

  xd    xd  y      y  2

2

d d

Thus, the wheel-rail normal contact force in the MKP-YAW contact model is deduced to be

E N 2 1 2



 2E  4 1 2



yl

xd  y   a  y 

yr

xd  y   a  y 

  

 wr  y  a 2  y    x  xd  y  

yr a 

yr

  xd    xd  y      y  2

a 2  y   wr  y 

yl

yr a 

2

d  d

(34)

dy

a 2     2

a  

yl

dxdy

a 2     2

a  

yl

2

  xd    xd  y      y  2

2

d  d

In the case of multi-point contact, the wheel-rail normal force is deduced to be N

 2E 4 1 2





 i 1 y Nc

y

a 2  y   wr  y 

li

ri

a  

 j 1 y lj  a  Nc

y

rj

dy

a 2     2   xd    xd  y      y  2

2

(35)

d  d

It can be found from the above descriptions that the MKP-YAW contact model has the same expression as the MKP model in calculating no-yaw angle contact as xd(y)=0. Moreover, the two-dimensional wheel-rail contact geometry calculation is used in the MKP-YAW contact model, which reduces the computational effort of the three-dimensional geometry calculation from the previous model [23]. Compared with the MKP model, the MKP-YAW method has little effect on the calculation efficiency, although this method considers the influence of the yaw angle. A detailed description concerning the calculation efficiency of the MKP contact model can be found in Refs. [23].

3.3. Verification of the MKP-YAW model The MKP-YAW contact model is validated by calculating the wheel-rail non-Hertzian contact between three types of wheels and rail with different lateral displacements and yaw angles. The calculation results obtained by KVM are taken as references. The validation conditions are proposed in Section 2.2, and the simulation conditions with a yaw angle of zero are also supplemented to reflect the influence of the yaw angle. 3.3.1 Verification of the contact patch According to Remark 2.3, the overlap ratio of the wheel-rail contact patch obtained by the MKP-YAW

contact model and KVM is applied to evaluate the computational accuracy of the shape of the contact patch.

(a)

(b)

(c) Figure 11. Overlap ratios for the wheel-rail contact patch obtained by MKP-YAW and KVM. (a) First type of wheel, (b) second type of wheel, and (c) third type of wheel.

Figure 11 shows the overlap ratios of the wheel-rail contact patch for the three types of wheel-rail contact conditions. It can be seen from Fig. 11(a) that the overlap ratio for calculating the standard wheel-rail contact is very high, and the minimum overlap ratio exceeds 92%. The overlap ratio for calculating the worn wheel and rail contact is commonly high but low for some particular calculation conditions. The MKP-YAW contact model uses the constant reduction factor  to calculate the contact area, which results in inaccurate estimation of the contact patch for certain contact conditions. A more detailed discussion of the influence of the reduction factor on the contact patch estimation can be found in references [22,27]. It is worth noting from Fig. 11 that the yaw angle changes the lateral displacement-overlap ratio distribution. However, the minimum overlap ratios under conditions are generally the same with or without the yaw angle. The level of calculation error when using the MKP-YAW contact model to estimate the contact area with the yaw angle is similar to the level when using the MKP model to estimate the non-yaw angle contact area. 3.2.2 Verification of the wheel-rail normal contact force The wheel-rail normal contact forces calculated by the MKP-YAW contact model and KVM are compared in Fig. 12. The distributions of the normal forces induced by different lateral displacements and yaw angles obtained by the two wheel-rail contact methods are very similar. Similar to the characteristics of the contact patch shown in Fig. 11, the yaw angle has a significant influence on the lateral displacement-normal force distribution. All the normal force nephograms have a similar variation characteristics, showing that the lateral displacement—normal force function shifts towards the direction of large lateral displacement as the yaw angle increases.

(a1)

(a2)

(b1)

(b2)

(c1)

(c2)

Figure 12. Comparison of the wheel-rail normal contact forces obtained by the MKP-YAW contact model (left) and KVM (right). (a) First type of wheel, (b) second type of wheel, and (c) third type of wheel.

It can be seen from Fig. 12 that the distributions of the normal force are quite different for the contact conditions of different wheel types. The normal force slightly changes with the lateral displacement and the yaw angle for the standard wheel, while there is an obvious increased wheel-rail normal force area for the worn wheel. The reason for this difference is that there is conformal contact at the increased normal force area for the second wheel type case, and there is two-point contact for the third type case. Both the conformal contact and the two-point contact result in an increase in the wheel-rail contact patch and the contact stiffness and further increase the wheel-rail normal force. A more detailed description can be found in Ref. [23,25].

(a)

(b)

(c) Figure 13. The relative errors of the wheel-rail normal contact force obtained by the MKP-YAW contact model. (a) First type of wheel, (b) second type of wheel, and (c) third type of wheel.

Taking the wheel-rail normal force calculated by KVM as the reference, the relative errors of the results obtained by the MKP-YAW contact model are shown in Fig. 13. As seen from Fig. 13, the calculation errors of the MKP-YAW contact model are quite small for the contact conditions of the first type of wheel and the third type of wheel, less than 2% and 4%, respectively. The calculation error for the contact condition of the second type of wheel is relatively large in some cases, close to 20%. The poor estimation of the contact area is the main reason for the larger calculation error of the wheel-rail force for the second contact case. The calculation error of MKP-YAW, in most cases, is relatively small, and the abrupt change in the normal force that commonly appears in other non-Hertzian contact models [23] does not exist in this contact model. MKP-YAW, therefore, can exhibit excellent calculation accuracy and computational stability. Similar to the overlap ratio of the contact patch shown in Fig. 11, the yaw angle can change the lateral displacement-normal force distribution. However, the error of the wheel-rail normal force calculated by the MKP, without considering the yaw angle, is similar to that calculated by the MKP-YAW model considering the yaw angle. 3.3.3 Verification of the wheel-rail normal contact pressure Wheel-rail normal contact pressure is a key factor affecting wheel-rail wear and wheel-rail rolling contact fatigue (RCF) [35,36]. To evaluate the computational accuracy of the contact pressure in the MKP-YAW contact model, the relative errors of the maximum contact pressure and the overlap ratios of the wheel-rail contact pressure obtained by the MKP-YAW contact model and KVM are adopted and shown in Fig. 14.

(a1)

(a2)

(b1)

(b2)

(c1)

(c2)

Figure 14. The relative errors of the maximum wheel-rail normal pressure (left) and the overlap ratios of wheel-rail contact pressure (right) obtained by the MKP-YAW contact model and KVM. (a) First type of wheel, (b) second type of wheel, and (c) third type of wheel.

It can be seen from Fig. 14 that the relative errors of the maximum wheel-rail normal pressure and the overlap ratios show the opposite distribution, which means that the contact conditions with a larger relative error for the maximum wheel-rail normal pressure are consistent with those contact conditions with lower overlap ratios of normal pressure. The relative errors of the maximum wheel-rail normal pressure are generally larger than the errors of the wheel-rail normal force that are shown in Fig. 13. The maximum relative error of the maximum wheel-rail normal pressure exceeds 20% for the contact conditions of the first and third types of wheels and even exceeds 40% for the contact conditions of the second type of wheel. The overlap ratios of the wheel-rail contact pressure shown in Fig. 14 are similar to the overlap ratios for the wheel-rail contact patch that are shown in Fig. 11. The wheel-rail contacts of the standard wheel (the first type of wheel) achieve high overlap ratios, while the overlap ratios of the certain wheel-rail contact conditions of the worn wheels are slightly lower. The minimum overlap ratio of the normal pressure occurs at the wheel-rail contact of the second type of wheel with the minimum overlap ratio of the contact patch and the maximum calculation error of the normal force and the maximum normal pressure. Similar to the contact area and the normal contact force, the computational accuracy level of the normal contact pressure does not decrease with increasing yaw angle in the MKP-YAW contact model. The above analyses indicate that by using the approximate expression (Eq. (15)) based on the basic assumptions, the MKP-YAW contact model can maintain the computational accuracy at the high level of the MKP model when calculating the non-yaw angle contact in the simulation of the wheel-rail non-Hertzian contact with yaw angle considered.

4. Conclusions A simplified model is presented for solving the wheel-rail normal contact problem under the influence of the yaw angle. Three basic assumptions relating to the approximate expression of the wheel-rail contact area and the normal contact pressure considering the yaw angle are first proposed. The computational accuracy of the approximate expression is validated by comparison with the contact behaviour obtained by KVM. The contact behaviours between the standard Chinese rail profile CHN 60 and three types of wheel profiles with different wheelset displacements and yaw angles are adopted in the validation. The calculation results indicate that the basic assumptions can be used to describe the wheel-rail contact patch and normal pressure with relatively high computational accuracy, and the approximate expression based on the basic assumptions is suitable for the simplified wheel-rail normal contact analysis. The validated approximate expression is applied to the MKP contact model, which is an approximate wheel-rail non-Hertzian normal contact model, to obtain the improved contact model—the MKP-YAW contact model. Compared with the MKP contact model, the MKP-YAW model can also consider the influence of the yaw angle. The computational accuracy, in terms of the wheel-rail normal contact force, the shape of the contact patch, and the normal contact pressure calculated by the MKP-YAW model, is comprehensively evaluated by comparison with the results obtained by KVM. The estimation results show that the computational accuracy of our proposed MKP-YAW model in calculating the wheel-rail contact considering the yaw angle is similar to that of the MKP contact model in calculating the wheel-rail contact without considering the yaw angle. The MKP-YAW achieves good agreement with KVM in the calculation of the wheel-rail normal force, while having a relatively worse estimation of the wheel-rail contact shape and normal pressure. The experimental results from the test rigs can more directly reflect the accuracy of the approximate expression and the MKP-YAW model. A well-designed experiment on test rigs will be very helpful for verification. The validation of the proposed model on test rigs requires further study. In this paper, the MKP-YAW contact model can be regarded as an application example of the approximate expression in solving the problem of wheel-rail non-Hertzian contact with yaw angle. The approximate expression can also be applied to other approximate wheel-rail non-Hertzian normal contact models in future research.

Acknowledgments This work was supported by the Open Research Fund of State Key Laboratory of Traction Power of Southwest Jiaotong University [Grant No. TPL2011], National Natural Science Foundation of China [grant Nos. 51735012 and 11790283], and Key Research & Development Program of Jiangsu Province [grant No. BE2019613].

CRediT authorship contribution statement Yu Sun: Methodology, Software, Validation, Writing - original draft, Funding acquisition. Wanming Zhai: Conceptualization, Funding acquisition, Writing - review & Editing. Yunguang Ye: Formal analysis, Writing - Review & Editing. Liming Zhu: Funding acquisition, Visualization. Yu Guo: Project administration, Investigation, Writing - Review & Editing.

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References [1]

Escalona JL, Aceituno JF. Multibody simulation of railway vehicles with contact lookup tables. Int J Mech Sci 2019;155:571-582.

[2]

Zhai W, Han Z, Chen Z, Ling L, Zhu S. Train–track–bridge dynamic interaction: a state-of-the-art review. Vehicle Syst Dyn 2019;57(7):984-1027.

[3]

An B, Ma D, Wang P, Zhou J, Chen R, Xu J, Cui D. Assessing the fast non-Hertzian methods based on the simulation of wheel–rail rolling contact and wear distribution. Proc Inst Mech Eng Part F 2019. https://doi.org/10.1177/0954409719848592

[4]

Hertz H. Über die berührung fester elastische Körper und über die Harte. fur diereine und Angew Math 1882;92:156–71.

[5]

Johnson KL. Contact Mechanics. Cambridge: Cambridge University Press; 1985.

[6]

Meymand SZ, Keylin A, Ahmadian MA. survey of wheel–rail contact models for rail vehicles. Vehicle Syst Dyn 2016;54(3): 386-428.

[7]

Wriggers P. Finite element algorithms for contact problems. Arch Comput Methods Eng 1995;2:1–49.

[8]

Zhao X, Li Z. The solution of frictional wheel–rail rolling contact with a 3D transient finite element model: Validation and error analysis. Wear 2011;271(1):444-452.

[9]

Yang Z, Boogaard A, Wei Z, Liu J, Dollevoet R, Li Z. Numerical study of wheel-rail impact contact solutions at an insulated rail joint. Int J Mech Sci 2018;138:310-322.

[10] Vo KD, Tieu AK, Zhu HT, Kosasih PB. A 3D dynamic model to investigate wheel–rail contact under high and low adhesion. Int J Mech Sci 2014;85:63-75. [11] Pun C L, Kan Q, Mutton P J, et al. An efficient computational approach to evaluate the ratcheting performance of rail steels under cyclic rolling contact in service. Int J Mech Sci 2015;101:214-226. [12] Xu Y, Jackson R L. Boundary element method (BEM) applied to the rough surface contact vs. BEM in computational mechanics. Friction 2019;7(4):359-371. [13] Kalker JJ. Three-dimensional elastic bodies in rolling contact. Dordrecht, The Netherlands: Kluwer Academic Publishers; 1990. [14] Vollebregt E. A new solver for the elastic normal contact problem using conjugate gradients, deflation, and an FFT-based preconditioner. J Comput Phys 2014;257:333–351. [15] Vollebregt E, Segal G. Solving conformal wheel–rail rolling contact problems Vehicle Syst Dyn 2014;52(sup1): 455-468. [16] Burgelman N, Li Z, Dollevoet R. A new rolling contact method applied to conformal contact and the train–turnout interaction. Wear 2014;321:94-105. [17] Deng X, Qian Z, Li Z, et al. Applicability of half-space-based methods to non-conforming elastic normal contact problems. Int J Mech Sci 2017;126:229-234. [18] Paggi M, Bemporad A, Reinoso J. Computational Methods for Contact Problems with Roughness//Modeling and Simulation of Tribological Problems in Technology. Cham: Springer; 2020:131-178. [19] Linder Ch. Verschleiß von Eisenbahnrädern mit Unrundheiten [Dissertation Nr. 12342]. Zurich: ETH; 1997.

[20] Piotrowski J, Kik W. A simplified model of wheel/rail contact mechanics for non-Hertzian problems and its application in rail vehicle dynamic simulations. Vehicle Syst Dyn 2008;46(1–2):27–48. [21] Piotrowski J, Chollet H. Wheel–rail contact models for vehicle system dynamics including multi-point contact. Vehicle Syst Dyn 2005;43(6-7):455-483. [22] Liu B, Bruni S, Vollebregt E. A non-Hertzian method for solving wheel–rail normal contact problem taking into account the effect of yaw. Vehicle Syst Dyn 2016;54(9):1226-1246. [23] Sun Y, Zhai W, Guo Y. A robust non-Hertzian contact method for wheel-rail normal contact analysis, Vehicle Syst Dyn 2018;56(12):1899-1921. [24] Sun Y, Guo Y, Zhai W. Prediction of rail non-uniform wear–Influence of track random irregularity. Wear 2019;420-421:235-244. [25] Sun Y, Guo Y, Lv K, Chen M, Zhai W. Effect of hollow-worn wheels on the evolution of rail wear. Wear, 2019;436-437:203032. [26] Ayasse JB, Chollet H. Determination of the wheel rail contact patch in semi-Hertzian conditions. Vehicle Syst Dyn 2005;43(3):161-172. [27] Sh. Sichani M, Enblom R, Berg M. A novel method to model wheel–rail normal contact in vehicle dynamics simulation. Vehicle Syst Dyn 2014;52(12):1752-1764. [28] Vo KD, Zhu HT, Tieu AK, Kosasih PB. FE method to predict damage formation on curved track for various worn status of wheel/rail profiles. Wear 2015;322:61-75. [29] An B, Wen J, Wang P, Wang P, Chen R, Xu J. Numerical investigation into the effect of geometric gap idealisation on wheel-rail rolling contact in presence of yaw angle. Math Probl Eng 2019. https://doi.org/10.1155/2019/9895267 [30] Zhai W, Liu P, Lin J, Wang K. Experimental investigation on vibration behaviour of a CRH train at speed of 350 km/h. Int J Rail Transp 2015;3(1):1-16. [31] Burgelman N, Li Z, Dollevoet R. Effect of the longitudinal contact location on vehicle dynamics simulation. Math Probl Eng 2016. http://dx.doi.org/10.1155/2016/1901089 [32] Zhai W. Vehicle–Track Coupled Dynamics Theory and Applications. Singapore: Springer; 2019. [33] Xu L, Zhai W. Stochastic analysis model for vehicle-track coupled systems subject to earthquakes and track random irregularities. J Sound Vib 2017;407:209-225. [34] Chen G, Zhai WM. A new wheel/rail spatially dynamic coupling model and its verification. Vehicle Syst Dyn 2004;41(4):301-322. [35] Meghoe A, Loendersloot R, Tinga T. Rail wear and remaining life prediction using meta-models. Int J Rail Transp 2019; 1-26. https://doi.org/10.1080/23248378.2019.1621780 [36] Enblom R. Deterioration mechanisms in the wheel–rail interface with focus on wear prediction: a literature review. Vehicle Syst Dyn 2009;47(6):661-700.

Appendix. Calculation of correction coefficients  and  The correction coefficients of  and can be expressed as    W / L     L / W

with,

(A1)

 L  2a  0   W  yr  yl

(A2)

  0 , if W / L   1 1 /  0 , if W / L   1

 

(A3)

0  0.5837   A / B   0.1053   A / B   0.5184  A / B  2





W / L 2 , if W / L   1  A/ B   2  L / W  , if  L / W   1

1

(A4)

(A5)

GA