Effect of the first two wheelset bending modes on wheel-rail contact behavior

Abstract

The objective of this paper is to develop a new wheel-rail contact model, which is suitable for considering the effect of wheelset bending deformation on wheel-rail contact behavior at high speeds. Dummies of the two half rigid wheelset are introduced to describe the spacial positions of the wheels of the deformed wheelset. In modeling the flexible wheelset, the first two wheelset bending modes are considered. Based on the modal synthesis method, these mode values of the wheelset axle are used to solve the motion equations of the flexible wheelset axle modeled as an Euler-Bernoulli beam. The wheel is assumed to be rigid and always perpendicular to the deformed axle at the wheel centre. In the vehicle model, two bogies and one car body are modeled as lumped masses. Spring-damper elements are adopted to model the primary and secondary suspension systems. The ballasted track is modeled as a triple layer discrete elastic supported model. Two high-speed vehicle-track models, one considering rigid wheelset models and the other considering flexible wheelset models, are used to analyze the differences of the numerical results of the two models in both frequency and time domains. In the simulation, a random high-speed railway track irregularity is used as wheel-rail excitations. Wheel-rail forces are calculated and analyzed in the time and frequency domains. The results clarify that this new contact model can characterize very well the influence of the first two bending modes of the wheelset on contact behavior.

概要

研究目的

扩展动力学模型的分析频域, 建立能考虑轮对柔性的车辆轨道耦合动力学系统模型, 为研究轮轨磨耗的形成和发展以及轮轨噪声的来源提供基础。

创新要点

利用欧拉梁横向弯曲模型, 建立轮轴在垂直于轨道平面和平行于轨道平面内的弯曲振动模型; 建立考虑轮轴弯曲的轮对模型与轮轨接触模块之间的耦合关系, 进而研究轮轨接触行为受轮轴弯曲变形的影响。 研究方法: 1. 把轮轴模拟为欧拉梁, 左右车轮模拟为固结于轮轴上的质量块; 2. 假设左右车轮始终垂直于轮轴, 引入虚拟的两个半边刚性轮对模型, 建立轮轨接触模型和柔性轮对耦合的关系; 3. 基于多刚体车辆-轨道耦合动力学模型, 利用以上柔性轮对模型和此耦合关系, 建立考虑轮轴柔性的车辆-轨道耦合动力学模型。

重要结论

1. 建立的刚柔耦合的车辆-轨道耦合动力学模型能够有效地描述轮轴弯曲对轮轨接触行为的影响; 2. 在0–150 Hz的随机不平顺激励下, 多刚体模型和考虑轮对柔性的模型受到的轮轨力和轮轨接触点轨迹不同; 这主要是由第1阶弯曲模态被激发导致。

1 Introduction

High-speed railways are currently popular globally. However, there are some problems including passenger riding comfort, noise pollution, and even operational safety (Jin et al., 2013). Rail corrugation, rail welding irregularity, wheel burning, and wheel out-of-roundness (OOR) generate high-frequency components of the dy-namic wheel-rail contact forces that contribute signifi-cantly to the total wheel-rail contact forces (Nielsen et al., 2003), and reduce the life of the components of track and vehicle, such as wheels, rails, and fasteners. Rail grinding and wheel re-profiling are the most common measures that have been proved to be effective in con-trolling rail irregularities and wheel OOR. However, these measures lead to notably high maintenance costs. A lot of measurements at the sites and coupling vehicle-track dynamics modeling have been carried out to in-vestigate the mechanism and development of these phenomena. In the vehicle-track dynamics modeling, a rigid multi-body system is often adopted to simulate railway vehicles, based on several commercial codes available for the low-frequency domain, such as GENSYS, NUCARS, SIMPACK, and VAMPIRE. These computer programs are generally used to analyze railway vehicle dynamics responses at frequen-cies below 20 Hz, where the influence of rigid motions of the vehicle on wheel-rail contact forces is dominant (Nielsen et al., 2005). To analyze the vehicle dynamic responses at mid- and high-frequencies, the vehicle structural flexibility should be taken into account in the modeling. It is obvious that wheelset structural flexibil-ity has an influence on wheel-rail contact behaviors at mid- and high-frequencies. Different flexible wheelset models have been set up due to various motivations in the past (Chaar, 2007).

The methods applied to modeling flexible wheelset can be summarized as three major categories (Chaar, 2007). The first is a lumped model developed in a simple and convenient way, in which a wheelset is divided into several parts interconnected with springs and dampers. This model can describe the bending and torsional mo-tions of the wheelset with only a few degrees of free-dom, which could not be applied to studying wear phenomena on wheel treads or rails (Popp et al., 1999). The second is a continuous model developed by Szolc (1998a; 1998b), in which the wheelset axle was modeled as a beam, and two wheels and brake discs were mod-eled as rigid rings attached to the axle through a mass-less, elastically isotropic membrane. The model can characterize the wheelset dynamic behavior in the fre-quency range of 30–300 Hz. In the model proposed by Popp et al. (2003), the wheelset axle was considered as a 1D continuum, having the properties of a bar, a torsional rod, and a Rayleigh beam. The wheel was considered as a 2D continuum, having the proper-ties of a disc and a Kirchhoff plate. The third was de-veloped based on finite element method (FEM), which simulates wheelset flexibility more realistically than the first two categories of model. The wheelset modes and corresponding natural frequencies were obtained through the modal analysis of the finite element (FE) model by using the commercial software, and they were input into the simulation by means of the commercial codes (SIMPACK, NEWEUL) (Meinders and Meinker, 2003) or some non-commercial multi-body dynamic system codes. The non-commercial code developed by Fayos et al. (2007) and Baeza et al. (2008; 2011) introduced the Eulerian coordinate system to replace the Lagrangian coordinate system in the flexible wheelset modeling. In this way it is convenient to obtain the motion of fixed physical nodes, and consider the inertial effect due to wheelset rotation. Relying on cur-rent computing power, it is feasible to use FEM to con-sider the effect of flexible wheelset in modeling a rail-way vehicle coupling with a track.

Regarding the wheel-rail contact treatment in con-sidering flexible wheelset influence, wheel-rail rolling contact condition is simplified based on different prior assumptions, especially in the detection of wheel-rail contact points. This is the prerequisite for the calculation of wheel-rail creepages and contact forces. Baeza et al. (2011) neglected the effect of the high-frequency de-formation and the deviation of a rotating flexible wheelset rolling over a flexible track model on the wheel-rail contact point in the investigation into the effect of the rotating flexible wheelset on rail corruga-tion. Through the detailed calculation Kaiser and Popp (2006) found that the contact point was in the location where the wheel and the rail had positive penetration maxima, and the penetration direction was orthogonal to the common tangent plane of the wheel and the rail before their deformations. A linear wheel-rail contact model was proposed and used to carry out the detection of wheel-rail contact point and the contact zone’s normal direction (Andersson and Abrahamsson, 2002). In the detection, the functions were created using a first-order Taylor expansion around a reference state described by a group of parameters which represent a configuration, in which the train was in static equilibrium and the wheel and the track were free from geometric imperfections. The advantage of this approach is that the contact position and orientation in each time step can be calculated by interpolation replacing iterations, which results in a low computa-tional cost. But the approach is only suitable for the case that the effect of all the parameters is very small on the contact point position and the contact patch orientation around the references is in static equilibrium. The wheel-rail contact point position and the contact patch orientation greatly depend on parameters, such as the curvatures of wheel and rail. In (Torstensson et al., 2012; Torstensson and Nielsen, 2011), the contact point detection was done before the simulation and used in the subsequent time integration analysis in the form of look-up table. The commercial software GENSYS al-lows for such calculations using the pre-processor KPF (from Swedish contact point function). In the KPF, the location and orientation of the contact patch were as-sumed to be dependent only on the relative displacement in the lateral direction between the wheelset and the rails, and hence the influence of the wheelset yaw angle was not taken into account. In some other papers de-tailed discussions on the wheel-rail contact model were omitted. In this study, the wheel-rail contact model considering the effect of wheelset flexibility (Zhong et al., 2013; 2014) is further improved and the new contact model is suitable for the analysis on the effect of the local higher-frequency deformation of the wheels on the wheel-rail contact behavior.

2 Vehicle-track coupling dynamic system

A flexible wheelset model (to be illustrated in Sec-tion 2.1) and a suitable wheel-rail contact model (to be discussed in Section 2.2) are integrated into the vehi-cle-track coupling dynamic system model. All parts of the vehicle system, except for its four wheelsets, are considered as rigid bodies. The primary and secondary suspension systems of the vehicle are modeled with spring-damper elements. A triple layer model of discrete elastic support is adopted to simulate the ballasted track. The rails are modeled as Timoshenko beams. The sleepers are modeled as rigid bodies and the ballast model consists of discrete equivalent masses. The equivalent spring-damper elements are used as the con-nections between the rails and the sleepers, the sleepers and the equivalent ballast bodies, and the ballast bodies and the roadbed. Fig. 1 shows the vehicle-track coupling dynamic system model. The equations of motion of each component of the vehicle excluding wheelsets and the track are illustrated in detail in (Xiao et al., 2007; 2008; 2010). The parameters and their values describing the dynamic models are given in Appendix A.

2.1 Flexible wheelset model

The wheelset structural flexibility is considered by modeling the wheelset axle as an Euler-Bernoulli beam in two planes, one perpendicular to the track centerline and the other parallel to the track level. The crossing effect of the bending deformations in the two planes is ignored. In the first two bending modes obtained using the modal analysis of the FE model of a wheelset, two wheels have little deformation (Fig. 2), and their frequencies are in the available frequency range (0–500 Hz) of an Euler-Bernoulli beam model. Therefore, two wheels can be treated as rigid bodies in this study.

There are two force systems acting on the wheelset, one is the wheel-rail contact forces and the other is the forces of the primary suspension system (Fig. 3).

In Fig. 3, O _fL and O _fR are the left and right points on the axle, respectively, where the primary suspension force systems are applied. O _CL and O _CR are the left and right contact points of wheel-rail, respectively. O indicates the origin of the coordinate system O-XYZ that is a coordinate system with a translational motion along the tangent track centerline at operational speed. If the speed is constant, this coordinate system is an inertial coordinate system, and therefore regarded as an absolute coordinate system (geodetic coordinate system).

To analyze the axle’s deformation, the force systems from wheel-rail interaction acting on the left and right wheel treads are translated to the nominal circle centers O _L and O _R, respectively, and extra moments are produced in the procedure of translating contact forces. Thus, the force systems acting on the axle in the two planes are obtained in Fig. 4.

Fig. 4

Force analysis diagram in the plane O-YZ (a) and in the O-XY plane (b)

The notations of the variables and symbols are de-fined in Table 1. The subscript p denotes the primary suspension, the subscripts x, y, and z denote X-, Y-, and Z-direction, respectively, and A denotes the axle.

Table 1 The notations of the variables

The differential equation for the flexural vibration of an Euler-Bernoulli beam (the axle) in the plane O-YZ is written as

((1))

where

((2))

((3))

The force analysis diagram of the two wheels in-cluding the D’Alembert forces is shown in Fig. 5, based on which differential equations of motion of the two wheels are written as

((4))

((5))

Note that the lateral accelerations of the wheels are assumed to be the same as the wheelset axle so there is no relative motion between wheels and axle.

Substituting the expressions of F _A(L,R)z and M _A(L,R)x obtained through Eqs.(4) and (5) into Eqs. (2) and (3), respectively, we can obtain:

((6))

((7))

((8))

where

((9))

Consider a solution of Eq. (8) in the form:

((10))

Using the calculus of variation (Qiu et al., 2009), the modal function satisfies:

((11))

((12))

((13))

where δ _ij is the Kronecker delta. For i=j, Eq. (11) can be written as

((14))

To obtain the mode shape functions with the wheelset axle modeled as a uniform Euler-Bernoulli beam carrying two particles (wheels), the segment of the beam from the left end to the first particle is referred to as the first portion, in between the two particles as the second portion and from the second particle to the right end as the third portion. The beam mode shape will be the superposition of the mode shapes of the three por-tions. The derivation of the mode shape functions is presented in Appendix B. The first three modes have the frequencies of f ₁= 111Hz, f ₂=245Hz, and f ₃=547Hz, respectively. These mode shape functions are normalized so as to satisfy Eq. (14), as shown in Fig. 6. The third mode is not in the fre-quency range of 0–500 Hz where the Euler-Bernoulli beam is available to analyze the system. Hence, the effect of the first two modes on dynamic responses is conducted in this study.

Fig. 6

The first three bending mode shapes of the wheelset

According to the modal analysis, we let the solution of Eq. (8) have the form:

((15))

Substituting Eq. (15) into Eq. (8), the differential equation can be written as

((16))

Multiplying both sides of Eq. (16) by $U_{zj}$ and integrating over the domain 0<y<L, we can obtain:

((17))

Using the orthogonality of the modal shape function as expressed in Eqs. (11) and (13), Eq. (17) can be written as

((18))

where

((19))

Eq. (18) can be expressed as

((20))

Eq. (20) can be written in the matrix form:

((21))

where

((22))

The explicit integral method illustrated in (Zhai, 2007) is used to obtain the vector $\left[ {\ddot q_{zj} } \right]$ of each acceleration coordinate.

For the vibration in the plane YOX, the dif-ferential equation expressed with respect to $\left[ {\ddot q_{xy} } \right]$ can be written as

((23))

The derivation of Eq. (23) is similar to that of Eq. (20) and omitted here. Eq. (23) can be expressed in matrix form:

((24))

where

((25))

2.2 Wheel-rail contact model

As mentioned in Section 2.1, the main concern in this work is the wheelset axle bending. The wheels are assumed to be rigid and their nominal rolling circles are always perpendicular to the deformed wheelset axle at their interference fit surfaces. Fig. 7 shows that the flexible wheelset moves from its initial reference state (O ₁(t ₁)) to its t ₂ status (O ₂(t ₂)), which is described in the plane of O-YZ. O ₁ is the center of the un-deformed wheelset at t ₁, and O ₂ is the center of the deformed wheelset at any time t. O ₁ O ₂ is the displacement vector of the wheelset center due to its rigid motion, and ϕ _R1 is the roll angle due to the wheelset rigid motion. The auxiliary line, $A_{\text{L}}^0 A_{\text{R}}^0$, is the central line of the un-deformed wheelset axle, $A_{\text{L}}^1 A_{\text{R}}^{\text{1}}$ is obtained by moving $A_{\text{L}}^0 A_{\text{R}}^0$ from O ₁(t ₁) to O ₂(t ₂), and $A_{\text{L}}^2 A_{\text{R}}^2$ is obtained through rotating $A_{\text{L}}^1 A_{\text{R}}^1$ by ϕ _R1. $A_{\text{L}}^2 A_{\text{R}}^2$ is actually the central line of the rigid wheelset axle at t ₂. Fig. 7 shows that the wheels are assumed to be rigid and always perpendicular to the deformed axle line at their connections at any time t ₂.

Fig. 7

A flexible wheelset moving from its initial reference state (O ₁(t ₁)) to its any status (O ₂(t)) in the plane of O-YZ

To clearly describe the new wheel-rail contact model, the dummies of the two rigid half wheelsets, as shown in Fig. 8, are employed to describe wheel-rail rolling contact behavior affected by the wheelset bending. The two dummies are indicated by DWL and DWR, respectively, and the wheels of the DWL and the DWR are assumed to overlap the left and right wheels of the flexible wheelset, respectively, all the time, namely, the motion of the assumed rigid wheels of the flexible wheelset can be described by the DWL and the DWR (Fig. 8). ϕ _R2 is the roll angle of the right wheel due to the bending deformation of the flexible wheelset. It is exactly the included angle between the line $A_{\text{L}}^2 A_{\text{R}}^2$ and the axle line of the right wheel or the wheel of the DWR.

Fig. 8

The relationship between the two rigid half-wheelset dummies and the flexible wheelset

It is not difficult to calculate the wheel-rail contact geometry considering the effect of the flexible defor-mation of the wheelset or the local high-frequency deformations of the wheels if the spatial po-sitions of the DWL and the DWR are determined. De-termining the spatial positions of the DWL and the DWR involves calculating their motion parameters, such as the lateral displacements of the centers of the wheels of the DWL and the DWR, indicated by y _DWL and y _DWR, respectively, the vertical displacements, z _DWL and z _DWR, the roll angles, ϕ _DWL and ϕ _DWR, and the yaw angles, ψ _DWL and ψ _DWR. These parameters are key to calculating the contact geometry of the flexi-ble wheelset in rolling contact with a pair of rails by using this new wheel-rail contact model. This will now be demonstrated in detail.

Fig. 8 describes the motion of the DWL and the DWR influenced by the wheelset bending and its rigid motion in the O-YZ plane only. After the rigid wheelset moves with the center displacement of O ₁ O ₂ and the rolling angle of ϕ _R1 in the O-YZ plane of the global reference, O-XYZ, its center position O ₁(t ₁) reaches the position O ₂(t ₂) and $A_{\text{L}}^0 A_{\text{R}}^0$ reaches (or becomes) $A_{\text{L}}^2 A_{\text{R}}^2$. Note that the vector O ₁ O ₂ and the roll angle ϕ _R1 around axis X are described in the O-YZ plane. The dash-dot line $A_{\text{L}}^1 A_{\text{R}}^1$ is through point O ₂(t ₂) and parallel to $A_{\text{L}}^0 A_{\text{R}}^0$. From Fig. 6, it is obvious that the rolling angle of the DWR caused by the wheelset rigid motion is just ϕ _R1, and that caused by the wheelset bending deformation is ϕ _R2, so the total rolling angle of the DWR is ϕ _DWR ϕ _R1+ϕ _R2, as shown in Fig. 6.

In addition, the displacement of the DWR is the vector O ₁ O _3R, which could be written as

((26))

In Fig. 8, the vector O ₂ A ₁ is parallel to O _3R A ₂ with the same length l ₀. l ₀ is actually the distance between the center of the wheel nominal circle and the center of the un-deformed wheelset. The vector O ₂ O _3R is parallel to A ₁ A ₂, with the same length. Thus, O ₁ O _3R can be written as

((27))

Moreover, the vector O ₂ A ₁ is described by {x ₁ y ₁ z ₁}[i j k]^T in O-XYZ, and can be obtained by rotating the vector {0 l ₀ 0}[i j k]^T (coinciding with the line $A_{\text{L}}^1 A_{\text{R}}^{\text{1}}$) about the X-axis by ε _DWR. O ₂ A ₁ is written as

((28))

The curve $\widehat{B_{\text{L}} B_{\text{R}}}$ (Fig. 6) is the deformed axle center line of the wheelset, which does not consider the influence of the rotation caused by the wheelset rigid motion. The point B _R is the center of the right nominal circle. The axle center line $(\widehat{O_{\text{2}} A_{\text{2}} })$ of the deformed wheelset, can be obtained by rotating $\widehat{B_{\text{L}} B_{\text{R}}}$ about the X-axis by ϕ _R1. According to the definition of the curve $\widehat{B_{\text{L}} B_{\text{R}}}$, the vector O ₂ B _R is defined as

((29))

where {Δx ₂ Δy ₂ Δz ₂}[i j k]^T is the displacement vector of the center of the right nominal circle due to the axle bending. Then the vector O ₂ A ₂ is defined as {x ₃ y ₃ z ₃}[i j k]^T, and can be written as

((30))

which is obtained according to the relationship between $\widehat{O_{\text{2}} A_{\text{2}}}$ and $\widehat{O_{\text{2}} B_{\text{R}}}$, or $\widehat{O_{\text{2}} A_{\text{2}}}$ obtained by rotating $\widehat{O_{\text{2}} B_{\text{R}}}$ by ϕ _R1. The wheelset center displacement vector O ₁ O ₂ is defined as {x ₀ y ₀ z ₀}[i j k]^T.

Substituting Eqs. (28) and (30) and the expression of O ₁ O ₂ into Eq. (27), the vector O ₁ O _3R can be written as

((31))

Similarly, when considering the wheelset bending deformation in the plane O-XY, the vector >O ₁>O _3R should be given as

((32))

where ψ _R1 and ψ _R2 are the yaw angles caused by the rigid motion and the bending deformation in the plane O-XY, respectively.

ψ _DWR=ψ _R1+ ψ _R2 is the total yaw angle of the DWR. Similarly, the position of the DWL can be obtained. When the positions of the two dummies are known at t ₂, the wheel-rail contact geometry can be calculated. Then the positions of the wheel-rail contact points are easily found and the wheel-rail contact forces can be calculated. The normal wheel-rail contact forces are calculated by the Hertzian nonlinear contact spring model, and the tangent contact forces and spin moments are calculated by means of the model by Shen et al. (1983). Compared with the con-ventional wheel-rail contact model (Wang, 1984; Zhai, 2007), this new wheel-rail contact model can charac-terize the independent high-frequency deformations of the two wheels of the flexible wheelset more conveniently.

3 Results and discussion

When a vehicle is running on an ideal track, it is only excited by sleepers. Note that the “flexible” wheelset model used in this section denotes the model considering the first two bending modes. The dynamic system with flexible wheelset models is used in the simulation on an ideal track at the speed of 300km/h. Figs. 9a and 9b show the vertical forces in the frequency domain in steady and unsteady stages, respectively. In the unsteady stage, the peaks appear not only at a set of harmonic frequencies nf _s (n=1, 2, 3, …) produced by passing sleeper but also at f _b1, while the influence of the second bending mode is small since there is no peak at f _b2. In the steady stage, the contribu-tion of the component at f _b1 is weak-ened and only the peaks at nf _s (n=1, 2, 3, …) remain. These results are reasonable because when a system comes to a steady stage, its responses only contain the component at the excitation frequency.

Fig. 9

The vertical contact force in steady (a) and unsteady (b) stages

Based on a large range of site measurements, the components of roughness on rails mostly appear in the range of 1–20 m. The natural frequencies of the first two bending modes are below 250Hz, meaning the available frequency of this model is limited. Therefore, the components of the random irregularity on the rails are mainly in the frequency range of 0–150 Hz at the speed of 300 km/h. Fig. 10a presents the local section of 900–950m in the time domain, and Fig. 10b shows the irregularity in the frequency domain. Note that the results below are from the steady stage.

Fig. 10

Random irregularity in the time domain (a) and frequency domain (b)

Figs. 11 and 12 show the wheel-rail contact forces acting on the rigid and flexible wheelsets in the time and frequency domains, respectively. As shown in Fig. 11a, the average of the oscillation of the lateral contact force acting on the flexible model is a little smaller than that on the rigid wheelset model, and the shapes of the os-cillation are different. As shown in Fig. 11b, the vertical contact forces acting on the two models oscillate around a similar average, while their shapes are different. These differences are caused by the wheelset flexibility.

Fig. 11

Lateral contact forces (a) and vertical contact forces (b) in the time domain

Fig. 12

Lateral contact forces (a) and vertical contact forces (b) in the frequency domain

In the frequency domain, the distributions of the components contained in both the lateral and vertical contact forces are in the excitation frequency range of the random irregularity. A peak at frequency 2f _s appears in Figs. 11a and 11b. The contribution of the component at frequency f _s is overwhelmed by the effect of the irregularity. In addition, the uniform distribution in 0–150 Hz of the irregularity results in the non-uniform distribution of contact forces. As shown in Figs. 12a and 12b, the components in 80– 150 Hz are higher than those in 0–80 Hz. This shows that under this present irregularity, this dynamic system is more sensitive to the excitation in 80– 150Hz than to those in 0–80 Hz.

In the frequency domain, the component at f _b1 of the lateral contact force acting on the flexible model is a little larger than that on the rigid model, as marked using the arrow in Fig. 12a. This shows that the first bending mode is excited, and the availability of the model to characterize the wheelset bending is proved. However, there is no evident differ-ence at f _b1 for vertical contact forces acting on the two models. This shows that the wheelset bending deformation has a stronger effect on the lateral contact force than on the vertical contact force.

The wheel-rail contact force is affected by the posi-tion of the lateral contact points. Figs. 13a and 13b show the oscillations of the contact points in lateral direction described in the body coordinate system attached to the rail cross-section in the time and frequency domains, respectively. The average of the magnitudes of oscilla-tion of the contact points on the flexible model in the time domain is larger than that on the rigid model. This is caused by the wheelset bending. Moreover, it can weaken the relative movement between rail and wheel caused by the irregularity. Therefore, it is one cause of the smaller average of the lateral contact force acting on the flexible model (Fig. 11a). As shown in Fig. 13b, the difference of the components between the two models at f _b1 is evident. This explains the difference in the time domain (Fig. 13a) and again shows the effectiveness of the proposed model.

Fig. 13

Oscillation of contact points in lateral direction in the time domain (a) and frequency domain (b)

4 Conclusions

In this study a new wheel-rail contact model is integrated into the high-speed vehicle-track coupling dynamics system model, which takes into account the effect of wheelset structural flexibility. Based on the new vehicle-track model the effect of the first two bending modes of the wheelset on wheel-rail contact behavior is analyzed under the random irregularity in a frequency range of 0–150 Hz. The numerical results of the rigid wheelset model and the flexible wheelset model are compared in detail. The following conclusions can be drawn from the results:

• 1. The present vehicle-track model considering flexible wheelsets can very well characterize the effect of the flexible wheelset on wheel-rail dynamic behavior.

• 2. Under the excitation, the shapes of the oscillations of the wheel-rail contact forces and contact points for the new and conventional vehicle-track models are differ-ent. The difference is caused by the excited first bending mode of the wheelset.

For future work, the first improvement to be con-sidered is to model a wheelset using the FEM or the Timoshenko beam theory to broaden the model’s available frequency range. This could allow it to help investigate the mechanisms behind the generation and development of wheel-rail wear and noise.

Appendices

Appendix A

The vehicle notations and track parameters are given in Table A1.

Table A1 The vehicle notations and track parameters

Appendix B

The axle is modeled as a uniform Euler-Bernoulli beam carrying two particles (wheels). The segment of the beam from the left end to the first particle is referred to as the first portion, in between the two particles as the second portion and from the second particle to the right end as the third portion. The beam mode shape will be the superposition of the mode shapes of the three portions. The mode shape of each portion has four constants of integration, i.e., a total of 12 for the three portions. It is necessary to satisfy: the boundary conditions; continuity of deflection and continuity of slope at the two ‘locations’; and compatibility of bending moments and compatibility of forces acting on the two particles.

Here we take the calculation of the mode shape functions in the plane O-YZ as an ex-ample. Fig. B1 shows a uniform Euler-Bernoulli beam O ₁ O ₃ of flexural rigidity EI _x , and length (R ₁+R ₂+R ₃)L carrying the first particle of mass m _w at axial coordinate R ₁ L from O ₁ and the second particle of mass m _w at axial coordinate R ₃ L from O ₃.

Fig. B1

Coordinate systems attached to the three sections of the wheelset axle

To write the equations of transverse vibrations of the system, three coordinate systems are chosen with origin at O ₁, O ₂, and O ₃. The choice of these coordinate systems has some algebraic advantages. In the text, the subscripts k=1, 2, and 3 refer to the first portion, the second portion, and the third portion of the beam, respectively. For free vibration of the beam at frequency, if the amplitude of vibration of the beam is U _zk (y _k ) at axial coordinate >(y> _k (in the range 0<y _k <CR _k L), then based on the Euler-Bernoulli bending theory, the bend-ing moment M _xk (y _k ), the shearing force Q _zk (y _k ), and the mode shape differential equation for the three portions are

((B1))

To express these equations in dimensionless form, one defines the dimensionless axial coordinate Y _k , amplitude Z _k (Y _k ), operator Dⁿ, dimensionless bending moment M _xk (Y _k ), shearing force Q _zk (Y _k ), and a dimensionless natural frequency Ω as follows:

((B2))

Therefore, Eq. (B1) can be expressed in the dimensionless form:

((B30)

Consider the solution of the previous equation as

$$\eqalign{ & Z_k (Y_k ) = C_{k1} \sin (\alpha Y_k ) + C_{k2} \cos (\alpha Y_k ) \cr & \;\;\;\;\;\;\;\;\;\;\;\;\; + C_{k3} \sinh (\alpha Y_k ) + C_{k4} \cosh (\alpha Y_k ). \cr}$$

((B4))

There are 12 unknown constants C _ki (i=1, 2, 3, 4) for the three segments.

For free vibration the D’ force and moment acting on the left wheel is $m_{\text{w}} \omega ^2 U_{z1} (R_1 L)$ and $J_{\text{w}} \omega ^2 U^\prime_{z1} (R_1 L)$ respectively (Fig. B2). Continu-ity of deflection and continuity of slope at O _wL together with compatibility of bending moments and compatibility of forces acting on the left wheel results in

((B5))

Fig. B2

Wheel diagrams including D’Alembert forces

The D’ force and moment acting on the left wheel is $m_{\text{w}} \omega ^2 U_{z2} (R_2 L)$ and $J_{\text{w}} \omega ^2 U^\prime_{z2} (R_2 L)$, respectively. Continuity of deflection and continuity of slope at O _wR together with compati-bility of bending moments and compatibility of forces acting on the right wheel results in

((B6))

Note that Eq. (B6) takes into account the contra di-rections of the axial coordinates y ₂ and y ₃. Eqs. (B5) and (B6) in di-mensionless form are

((B7))

((B8))

For the free boundary condition at the left end of the first portion and the right end of the third portion, the coefficients of the dimensionless mode shape functions satisfy

$$C_{11} = C_{13} ,\;C_{12} = C_{14} ,\;C_{31} = C_{33} ,\;C_{32} = C_{34} .$$

((B9))

Then one can write the dimensionless mode shape functions of the first and third portions as

$$\eqalign{ & Z_k (Y_k ) = C_{k1} \left( {\sin (\alpha Y_k ) + \sinh (\alpha Y_k )} \right) \cr & \;\;\;\;\;\;\;\;\;\;\;\;\; + C_{k2} \left( {\cos (\alpha Y_k ) + \cosh (\alpha Y_k )} \right) \cr & \;\;\;\;\;\;\;\;\;\;\; = B_{k1} P_k (Y_k ) + B_{k2} V_k (Y_k ), \cr & \;\;\;\;\;\;\;\;k = 1,\;3. \cr}$$

((B10))

Substituting Eqs. (B10) and (B4) into Eqs. (B7), we can obtain:

((B11))

Then one can write the dimensionless mode shape function of the second portions as

$$Z_2 (Y_2 ) = B_{11} P_2 (Y_2 ) + B_{12} V_2 (Y_2 ),$$

((B12))

where

((B13))

The coefficients of sin(α Y ₂), cos(α Y ₂), sinh(α Y ₂), and cosh(α Y ₂) in Eq. (B4) (when k=2) correspond to those in the expression obtained by substituting Eq. (B13) into Eq. (B12), so we can obtain

$$C_{2i} = B_{11} P_{2i} + B_{12} V_{2i} ,\;\;\;\;i = 1,\;2,\;3,\;4.$$

((B14))

The coefficients of B ₁₁ and B ₁₂ in Eq. (B14) correspond to those in the expression by simplifying Eq. (B11), so we can obtain:

((B15))

((B160)

So far the mode shape functions of the second and third portions have four unknown constants in total. These four unknown constants can be calculated using Eq. (B8). The first two equations of Eq. (B8) can be written as

((B17))

One considers

((B18))

where

$$\eqalign{ & P_{31} = \frac{{{\text{D}}V_3 (R_3 )P_2 (R_2 ) + V_3 (R_3 ){\text{D}}P_2 (R_2 )}} {{{\text{D}}V_3 (R_3 )P_3 (R_3 ) - V_3 (R_3 ){\text{D}}P_3 (R_3 )}}, \cr & P_{32} = \frac{{{\text{D}}V_3 (R_3 )V_2 (R_2 ) + V_3 (R_3 ){\text{D}}V_2 (R_2 )}} {{{\text{D}}V_3 (R_3 )P_3 (R_3 ) - V_3 (R_3 ){\text{D}}P_3 (R_3 )}}, \cr & V_{31} = - \frac{{{\text{D}}P_3 (R_3 )P_2 (R_2 ) + P_3 (R_3 ){\text{D}}P_2 (R_2 )}} {{{\text{D}}V_3 (R_3 )P_3 (R_3 ) - V_3 (R_3 ){\text{D}}P_3 (R_3 )}}, \cr & V_{32} = - \frac{{{\text{D}}P_3 (R_3 )V_2 (R_2 ) + P_3 (R_3 ){\text{D}}V_2 (R_2 )}} {{{\text{D}}V_3 (R_3 )P_3 (R_3 ) - V_3 (R_3 ){\text{D}}P_3 (R_3 )}}. \cr} $$

((B19))

Using the last two equations of Eq. (B8), we can obtain:

((B20))

where

((B21))

Using the matrix form of Eq. (B20), one can obtain

((B22))

Eq. (B22) is the frequency equation, which is a transcendental equation. By using an iterative procedure based on linear interpolation, the first three natural fre-quencies are f ₁=111 Hz, f ₂=245 Hz, and f ₃=547 Hz, respectively.

The calculation of the coefficients of the three mode shape functions are demonstrated in detail in the fol-lowing.

The dimensionless mode shape functions can be written as

$$Z_k (Y_k ) = B_{11} P_k (Y_k ) + B_{12} V_k (Y_k ),\;\;\;k = 1,\;2,\;3.$$

((B23))

One may set the deflection of the first particle to be A and without loss of generality one may choose A=1, hence

$$A = 1 = Z_1 (R_1 ) = B_{11} P_1 (R_1 ) + B_{12} V_1 (R_1 ).$$

((B24))

From the above equation and Eq. (B20), one can obtain the following equations:

((B25))

Subsequently, substituting the last equations into Eq. (B23) (assuming k=1) one can obtain the dimensionless mode shape function of the first portion $Z_1 (Y_1 )\;(0 \leqslant Y_1 \leqslant R_1)$. Substituting Eq. (B25) into Eq. (B23) (assuming k=3) one can obtain the di-mensionless mode shape function of the third portion $Z_3 (Y_3 )\;(0 \leqslant Y_3 \leqslant R_3)$. By inserting Eqs. (B15) and (B16) into Eq. (B4) (assuming k=2) one can obtain the dimensionless mode shape function of the second portion $Z_2 (Y_2 )\;(0 \leqslant Y_2 \leqslant R_2)$.

Hence, the coefficient of the three mode shape functions can be calculated in Table B1.

Table B1 Coefficients of the three modes