Bifurcation analysis of a railway wheelset with nonlinear wheel–rail contact

Stability is a key factor for the operation safety of railway vehicles, while current work employs linearized and simplified wheel/rail contact to study the bifurcation mechanism and assess the stability. To study the stability and bifurcation characters under real nonlinear wheel/rail contact, a fully parameterized nonlinear railway vehicle wheelset model is built. In modeling, the geometry nonlinearities of wheel and rail profiles come from field measurements, including the rolling radius, contact angle, and curvatures, etc. Firstly, four flange force models and their effects on the stability bifurcations are compared. It shows that an exponent fitting is more proper than a quintic polynomial one to simulate the flange and works well without changing the Hopf bifurcation type. Then the effects of each term of the nonlinear geometry of wheel/rail contact on the Hopf bifurcation and limit circle bifurcation are discussed. Both the linear term and nonlinear term of rolling radius have a significant influence on Hopf bifurcation and limit point of circle (LPC) bifurcation. The linear critical speed (Hopf bifurcation point) and the nonlinear critical speed (LPC bifurcation point) changes times while within the calculated range of the linear term of the rolling radius. Its nonlinear term changes the bifurcation type and the nonlinear critical speed almost by half. The linear term of the contact angle and the radius of curvature of the wheel and rail profiles should be taken into consideration since they can change both the bifurcation point and type, while the cubic term can be ignored. Furtherly, the field measured wheel profiles for several running mileages are employed to examine the real geometry nonlinearities and the corresponding Hopf bifurcation behavior. The result shows that a larger suspension stiffness would increase the running stability under wheel wear.

Abstract Stability is a key factor for the operation safety of railway vehicles, while current work employs linearized and simplified wheel/rail contact to study the bifurcation mechanism and assess the stability. To study the stability and bifurcation characters under real nonlinear wheel/rail contact, a fully parameterized nonlinear railway vehicle wheelset model is built. In modeling, the geometry nonlinearities of wheel and rail profiles come from field measurements, including the rolling radius, contact angle, and curvatures, etc. Firstly, four flange force models and their effects on the stability bifurcations are compared. It shows that an exponent fitting is more proper than a quintic polynomial one to simulate the flange and works well without changing the Hopf bifurcation type. Then the effects of each term of the nonlinear geometry of wheel/rail contact on the Hopf bifurcation and limit circle bifurcation are discussed. Both the linear term and nonlinear term of rolling radius have a significant influence on Hopf bifurcation and limit point of circle (LPC) bifurcation. The linear critical speed (Hopf bifurcation point) and the nonlinear critical speed (LPC bifurcation point) changes times while within the calculated range of the linear term of the rolling radius. Its nonlinear term changes the bifurcation type and the nonlinear critical speed almost by half. The linear term of the contact angle and the radius of curvature of the wheel and rail profiles should be taken into consideration since they can change both the bifurcation point and type, while the cubic term can be ignored. Furtherly, the field measured wheel profiles for several running mileages are employed to examine the real geometry nonlinearities and the corresponding Hopf bifurcation behavior. The result shows that a larger suspension stiffness would increase the running stability under wheel wear.
Keywords Railroad vehicle dynamics Á Wheel/rail contact Á Bifurcation Á Hunting stability Á Flange force

Introduction
Since hunting stability determines the operation speed limit and the running safety of a railroad vehicle, it is an important topic in the theoretical study for product design. Hunting motion is a kind of self-excited, sustained lateral and yaw vibration of wheelsets and bogies, occasionally interconnected also with severe car body vibration. The lowest speed, at which the hunting appears, is referred to as critical speed. With a reasonable safety margin considered, the field operation speed of the train would be determined to ensure running safety and avoid alarm conditions. The hunting motion is monitored all the time as the train runs. Once its magnitude reaches the limit, an alarm would be triggered and measurement will be taken to avoid severe vibrations on the vehicle as well as large wheel/rail interacting forces, even introduces safety issues, i.e., track shift and derailment [1][2][3].
A free wheelset would be unstable and hunting motion is observed as long as the running speed is not zero. To increase hunting stability, suspensions are applied on its two ends to constrain it to the bogie frame and/or the car body. Then the bifurcation characteristics of hunting motion are closely related to the wheel/rail interaction and suspension parameters. The bifurcation determines the critical speed of the vehicle as well as the amplitude of hunting motion, which are related to derailment safety. Wheel/rail interaction contributes to the basic dynamics of railway vehicles, and the contact patch and forces are influenced by the geometrical parameters of wheel and rail profiles as well as the wheelset motions concerning the rail. Due to the complexity, difficulty, and computational cost in deriving the analytical form of the governing equation of motion for the complete vehicle system involving nonlinear wheel/rail interaction and its solving, a linearized and simplified wheel/rail interacting algorithm is widely employed in the mechanism study of the bifurcation characteristics of hunting, see [4][5][6][7][8][9] for detail. The simplification and linearization mainly concentrate on the wheel/rail contact geometry and forces calculation; for instance, the wheel is simplified as a cone with a nonlinear stiffness representing the flange contact by ignoring the geometrical nonlinearities [9][10][11][12][13][14][15]. This leads to unrealistic stability results compared to a nonlinear equivalent conicity and wheel/rail forces, as discussed in Refs. [16,17].
The wheel/rail interaction geometry is highly dependent on the profile shape of wheel and rail, and the wear on wheel and rail further introduces geometry variance [3,18,19]. As discussed by Polach in Refs. [19] and [20], although the equivalent conicity for the wheelset amplitude of 3 mm, which is usually considered in railway practice, is the same for two pairs of wheel and rail profiles, they possess significantly different contact states leading to opposite signs of the nonlinearity parameter of the contact geometry. This can be even extended to the curvatures, contact angles that both enter into the contact patch evaluation of wheel/rail interaction.
Current research dealing with nonlinear stability assessment of railway vehicles often investigates simplified models; either a single free/suspended wheelset (see Refs. [7,[9][10][11][12][13]) or a truck/half-vehicle model (see Refs. [4,14,21,22]), respectively. Researches in [5,6,8,[18][19][20][21]23] employ a 2-axle wagon model, and the modeling technique for lateral stability study was reviewed in Ref. [24]. Those models can also be classified by how many components or DOFs are involved, such as a free or constrained/suspended wheelset model (one body with 2 DOFs) [7][8][9], a rigid-steering truck model (one body with 2 DOFs, a pretty large steering stiffness between the wheelsets and the truck frame) [4] or a soft-steering truck model [12] (three bodies with 6 DOFs, a small steering stiffness between the wheelsets and the truck frame), a simplified two-truck vehicle model (seven bodies with 14 or 17 DOFs) [5,8]. Using these models, the Hopf bifurcation point is obtained through evaluating the eigenvalues of the Jacobian matrix concerning certain parameters, for instance, the vehicle speed, while the nonlinear critical speed has to be retrieved from the Poincaré map and the phase portrait through time-domain analysis. However, the bifurcation analysis is seldom carried out for models with a large number of DOF combing using a nonlinear wheel-rail contact model and worn profiles.
The nonlinear stability performance on the straight and curved tracks has been investigated, for instance, by Lee and Cheng [5,6], by Zeng [8], and by Zboinski and Dusza [23]. Lee and Cheng [5,6] examined the linear critical speed on both tangent and curved tracks using a simplified-and linearized-vehicle model, and the model is restricted to consider the motion of the vehicle in the horizontal plane. Besides a time-domain integrating algorithm, the Lyapunov indirect method was employed in their work to perform stability analysis. Ref. in [7] pointed out that the wheel conicity, the mass and inertia of the wheelset, and the steering stiffness are critical to the Hopf bifurcation characteristics. The aforementioned models and solving algorithms are generalized methods for studying the linear critical speed of hunting and the resultants are acceptable if the vehicle system were properly linearized. Investigations of the influence of vehicle/track elasticity on bifurcation and critical speed have been carried out by Kaiser and Popp [25], Zhai and Wang [26], and Ling et al. [27].
The stability assessment in the railway industry requires the possibility to apply the stability analysis on large, complex wheel/rail interactions under realistic operating conditions. These conditions include the profiles of wheel and rail, the contamination on the contact surface, and the track stochastic irregularities, etc. This inevitably results in a large number of computations, and the field measured profiles of wheel and rail and realistic friction do lead to divergence between the nonlinear critical speed (Saddle-node bifurcation) and linear critical speed (Hopf bifurcation) as assessed in Ref. [28].
In this investigation, an improved railroad vehicle wheelset model is presented to study the nonlinear bifurcation characteristics with nonlinear wheel/rail interaction, which takes the geometrical nonlinearities of wheel and rail profiles into account. The wheel/rail interaction employs a mathematical expression to calculate the size of the Hertzian contact patch, in which the rolling radius, contact angle, and radius of curvature of wheel and rail profiles are fitted as functions of wheelset lateral displacement. Then the tangential creep forces of wheel/rail contact are analytically expressed as the Hertzian contact patch geometry, and lead to a continuous and fast calculation compared to a look-up table interpolation. The organization of the paper is as follows. Section 2 introduces the vehicle modeling, including the geometrical nonlinearities of profiles and the normal and tangential contact forces real-time calculation, as well as the solving algorithm for the governing differential equations of motion. In Sect. 3, the flange force models are compared to show their effects on the bifurcation behavior of the vehicle system, and the limitations are pointed out in the linearized vehicle model. In Sects. 4 and 5, the Hopf bifurcation and LPC bifurcations are discussed, respectively, regarding the geometrical nonlinearities of wheel and rail profiles. The vehicle stability under wheel wear is discussed in Sect. 6. The last section summarizes and concludes the main results.

Vehicle modeling
A free wheelset is constrained by the rail through wheel/rail contact and the suspensions at the axle ends, as in Fig. 1. The suspension consists of springs in parallel with viscous dashpots to simulate the connection between the wheelset end and the bogie frame (or car body). In Fig. 1, red arrows refer to the suspension forces and gravity, green arrows refer to the nonlinear wheel/rail contact forces, a three combined one-way black arrow forms a coordinate system attached to a body, and a two-end black arrow illustrates the geometrical dimension of suspension or a body. The origin of the coordinate system on wheelset locates at its centre of gravity, which simplifies the governing differential equations of motion of the vehicle system. Two degrees of freedom are introduced for the mechanism study of hunting stability, lateral displacement, and yaw rotation [9,11,13].

Geometrical nonlinearities of wheel and rail
The wheel/rail contact forces are determined by the wheel/rail contact patch (point), which is highly dependent on the geometrical profiles of wheel tread and railhead. The geometrical parameters around the contact point consist of two radii of curvature on wheel and rail, as depicted in Fig. 2, and the contact angle on a wheel, as in Fig. 3. Compared to existing mechanism studies on the bifurcation analysis for the railway vehicle, which assuming the wheel tread as a cone and ignoring the curvature effects, this work takes those geometrical nonlinearities into account when evaluating the normal and tangential contact forces. Figure 2 illustrates the geometrical parameters around the wheel/rail contact point. For the wheel profile, the radius of curvature along the running direction is also named the rolling radius; meanwhile, the other one   Fig. 1 Wheelset and primary suspension forces and wheel/rail forces along with the lateral direction changes with the contact location on the wheel. For the rail profile, the radius of curvature along the running direction is assumed as infinity, while the other one along the lateral direction varies with the location of the contact point.
To take the nonlinear wheel profile rather than a cone one into account, a general technique is to fit the profile as a high-order polynomial. For instance, the rolling radius of wheel r r,l (subscript r and l refers to the right and left wheel, respectively) is expressed as a function of the wheelset lateral displacement y w , in which, a 1 , a 2 , and a 3 are polynomial coefficients through curve fitting for a practical wheel profile. This equation can be generalized written as r r,l = f(y w ) to further concern various curve fitting techniques. Similar to the variant in rolling radius, the contact angle shall also be expressed as follows: where c 0 , c 1 , c 2 , and c 3 are polynomial coefficients for fitting the contact angle. The lateral radius of curvature of wheel and rail, respectively, is furtherly fitted as polynomials, in which, e 0 , e 1 , e 2 , and e 3 are polynomial coefficients for a wheel, whereas f 0 , f 1 , f 2 , and f 3 for rail.

Nonlinear Hertzian normal contact
Follow the Hertzian contact theory [29], two coefficients are assumed constant in the neighborhood of the contact point O c , and linked to the main local curvatures, in which R w1 = r and R r1 = ?.
Conventionally, e a and e b are defined as semi-axis of the elliptical contact patch demonstrated at the right-bottom in Fig. 2, and their values can be    [30][31][32], a reliable expression of the e b /e a as a function of A/B [0, ?] is proposed in [33] with e b /e a & (A/B) 0.63 ; then the area of the contact ellipse pe a e b is analytically expressed as, in which E is Young's modulus and m is the Poisson's ratio, assuming the same material for the rail and the wheel; N represents normal contact force. In this expression, the term in square brackets contains the material and geometrical constants, and the second term contains only the normal contact force. This expression has the advantage of being continuous, simple, and fast to calculate compared to the look-up table. Thus, the contact patch size is fully presented as the geometrical parameters of wheel and rail profiles, and e a e b is furtherly used to evaluate the following tangential contact force.

Nonlinear tangential/creep forces
The creep forces are expressed as nonlinear functions of the creepages, determined by the vehicle speed and wheelset motions concerning the rail. Figure 3 illustrates the wheel/rail interaction geometry, in which three coordinates system are defined. o t x t y t z t represents the track trajectory coordinate system of wheelset and defines the position of wheelset on track and the track orientations; o w x w y w z w represents the coordinate system attached to the wheelset center but not rotate with its axle axis; o cl x cl y cl z cl and o cr x cr y cr z cr represent coordinates at the left and right contact point, respectively. Coordinate transformations are needed among deriving the creepages and the wheel/ rail forces.
In this simplified vehicle model, the vertical forces on each wheel are assumed to be the wheel load under a quasi-static force equilibrium assumption. It is furtherly assumed that the vertical force only comes from the decomposition of the normal force, then the normal force is retrieved as, where p 0 is the wheel load with a constant value p 0 = (m c ? 2m f ? 4m w )g/8, / w is the roll angle of a wheelset which is calculated from the rolling radius difference between the left and right wheel. The lateral decomposition of the normal force will enter into the equations of motion of wheelset, The longitudinal, lateral, and spin creepages on the left and right contact point concern the wheelset displacements, and velocities are summarized as: where n xl;r ; n yl;r and n sl;r represent the longitudinal, lateral, and spin creepages, respectively, for the left and right wheel, in which the creepages are defined at the contact coordinate system. By introducing the track vibrations into preceding equations, the track irregularities can be included [27]. The modified Kalker's linear creep model [34] is used to derive the creep forces and moments, F xcl;r ¼ Àf 11l;r n xl;r F ycl;r ¼ Àf 22l;r n yl;r À f 23l;r n sl;r M zcl;r ¼ f 23l;r n yl;r À f 33l;r n sl;r where f 11 = Ge a e b C 11 , f 22 = Ge a e b C 22 , f 23 = G(e a e b ) 3/ 2 C 23, and f 33 = G(e a e b ) 2 C 33 are creep coefficients and dependent on the geometry parameters at the contact point, and C 11 , C 22 , C 23 , and C 33 are Kalker coefficients determined by e a /e b through a look-up table interpolation. To save the interpolating procedure for fast calculation, a polynomial fit of coefficients C ij is proposed in [33], in which e b /e a is defined in the notations of Eq. (6). Since a frequency-domain analysis is usually needed, properly analytical expressions regarding the wheel/ rail interaction are crucial to the computation accuracy of vehicle responses as well as efficiency. Equation (10) is furtherly modified based on Shen-Hedrick-Elkins' nonlinear creep theory which introduces the tangential creep force saturation concerning the friction [35,36], where l is the friction coefficient at the contact surface Transform the creep forces in the contact coordinate system to the track coordinate system, As a resultant force of creep forces and normal force, the total wheel/rail contact forces exerted on the wheelset is concluded as, Since a continuous representation of the creep forces is required in the frequency-domain analysis, the force saturation is neglected, while this can be realized in a time-domain analysis.

Differential equations of motion
The governing differential equations of motion for the lateral displacement y w and yaw angle w w of the wheelset are given, respectively, by, in which F xp and F yp refer to the longitudinal and lateral suspension forces on each wheelset side, respectively, and are presented as, The governing differential equations of motion for the vehicle system can be further post-processed into a matrix form, where M, C, and K are, respectively, the mass, damping, and stiffness matrix, F s is the generalized force vector, and x is the generalized coordinates vector of the system and defined as, The value of each parameter listed in Table 1 is taken for simulation. The value of each nonlinear term of wheel and rail is fitted from the contact calculation of a measured worn LMB10 wheel profile and CN60 rail profile used in China Railway.

Solving algorithm
To observe the bifurcation characteristics of this fully nonlinear wheelset model, both time-and frequencydomain analyses are performed. In the time-domain analyses, the geometrical nonlinearities of profiles and nonlinear creep forces can be completely concerned, whereas properly linearization and analytical expressions are needed in the pre-processing stage in the frequency-domain analysis. Equation (20) can be reexpressed as the following system of first-order differential equations, where x(t) is the vector of state variables. As well discussed, the equilibrium of railway wheelset on tangent tracks is x = 0. When it loses stability with the changing parameters, a Hopf bifurcation can be observed. The vehicle stability can be evaluated by the Lyapunov indirect method [5,8].
For any given parameter, such as the running speed v, one defines a determinant matrix A, where x 0 is the equilibrium state. The stability of the system can be assessed by the eigenvalues of matrix A. Besides, the first Lyapunov coefficient (FLC) can be calculated using the central manifold theorem and normal form to tell the Hopf bifurcation type, subcritical or supercritical.
The piecewise linear model is close to reality and widely used in ODE integral (time-domain integration of ordinary differential equations of motion); however, it is not suitable in bifurcation analysis since it is not continuous. The quintic polynomial model is continuous and suitable for symbolic analysis. However, simulation results in Table 2 show that it changes FLC so significantly that a subcritical Hopf bifurcation is observed instead of a supercritical bifurcation even though the Hopf bifurcation velocity (the linear critical speed, v cr_lin ) remains almost the same with or without the polynomial flange force considered. It is recommended to ignore the flange force if only the Hopf bifurcation is discussed. To calculate the limit circle (LC), a flange model is still needed. Thus, a new flange model should be taken instead of the polynomial flange force. An exponent model and reciprocal model are assumed as in Eqs. (26) and (27), respectively, F t ¼ p 1 e q 1 y w À p 1 e Àq 1 y w ð26Þ where the coefficients can be fitted through a MATLAB package called Curve Fitting Tool (cftool), p 1 = 0.0003906, q 1 = -1645, p 2 = -1, q 2 = 0.0112.
The flange forces calculated through each flange model are illustrated in Fig. 4. The exponent model is coincident with the piecewise model when the lateral displacement is within ± 11 mm and increases and decreases sharply once the lateral displacement gets out of the range. Test results in [13] show that the amplitude of the wheelset hunting motion is about 10 mm; thus the exponent model meets the requirement of simulation. The polynomial model shows a noticeable difference with others when the lateral displacement is within 9-11 mm. Within this range, the flange force has the same direction as the lateral displacement, i.e., it would contribute to the hunting motion, which is not following reality. The reciprocal model is coincident to the piecewise model when the lateral displacement is within ± 11 mm, but the flange force is much smaller than others with a lateral displacement of 11 mm; thus its applicability needs to be discussed furtherly.
A few terms are neglected in this section, such as the damping of primary suspension, the nonlinear wheel/rail contact parameters. Each case represents one type of Hopf bifurcation, namely supercritical (H-, with a negative FLC) or subcritical (H?, with a positive FLC). Calculated results are listed in Table 2. The reciprocal model gets a bit larger v cr_lin and smaller FLC in both cases. The exponent results are the same as the original model without the flange forces.
The polynomial flange model does not change v cr_lin , but it changes the bifurcation type. This is caused by the extremely large cubic term, a smaller FLC can be obtained with a smaller d 1 . Moreover, the cubic term of the polynomial flange model contributes to the instability of the system a lot, which is inconsistent with reality. ODE integral in Fig. 5 shows that the system experiences heavy hunting motion with an amplitude of 10 mm with a polynomial flange model, while the system return to its equilibrium at 0 without the flange forces.
The bifurcation map is illustrated in Fig. 6. For case 1, the magnitude of the limit cycle (y w ) of the original model without flange increases with velocity and reaches up to 13 mm, while that of the exponential model is significantly smaller. A period-doubling bifurcation of circles is observed at v = 57.5 m/s with reciprocal flange force, which implies that the reciprocal model is unappropriated.  From the bifurcation map, one can tell that the exponential flange model is much closer to reality than the others, the magnitude of LCs is much smaller than that without the flange model and the LC stays stable since the effect of the flange. In the following analysis, the exponential flange model is taken.

Analysis of Hopf bifurcation
With the default parameters in Table 1, a Hopf bifurcation occurs at point H: v = 51.70 with a FLC of -28.01, which has a little difference from the results from Sect. 3 since more nonlinear characters are included. To discuss the influence of the nonlinear wheel/rail contact geometry on Hopf bifurcation, a Hopf bifurcation analysis is conducted taking running speed and one of the terms of polynomials. The Hopf bifurcation curve is illustrated in Fig. 7, where Hstands for negative FLC, H? for positive FLC, GH for generalized Hopf bifurcation point with a FLC of zero. From Fig. 7a, one can tell, a 1 , the linear term of the rolling radius significantly affects the Hopf bifurcation speed, i.e., the v cr_lin . It decreases with a ratio inversely proportional to a 1 . Two generalized Hopf bifurcation points are detected at point GH: a 1 = 11.14 and a 1 = 38.14, which are clearly out of the reasonable value range. From Fig. 7b, one can tell that the v cr_lin varies little when a 3 changes a lot. A generalized Hopf bifurcation point is observed at point GH: a 3 = 102.47 with a second Lyapunov coefficient (SLC) of -30.01. A codimension two bifurcation is expected and the dynamic behavior is much more complex at this point [9]. Supercritical Hopf bifurcation occurs when a 3 \ 102.47, subcritical Hopf bifurcation occurs when a 3 [ 102.47. However, it is unclear whether the value of a 3 can be taken in reality.
As in Fig. 7c, the v cr_lin varies noticeably with c 1 . Two generalized Hopf bifurcation points are observed at point GH: c 1 = 164.23, c 1 = 1137.95 with SLC of 82.46 and 3.47 9 10 10 respectively. The latter is clearly out of a reasonable value range. The v cr_lin stays the same when c 3 changes a lot as in Fig. 7d. No GH point is observed within the calculated range.
As in Fig. 7e, The v cr_lin varies slowly with e 1 . Two generalized Hopf bifurcation points are observed at point GH: e 1 = 21,023.16, e 1 = -2152.34 with nonzero SLCs. The v cr_lin stays the same when e 3 changes a lot as in Fig. 7f. No GH point is observed within the calculated range. The v cr_lin varies slowly with f 1 in Fig. 7g. Two generalized Hopf bifurcation points are observed at point GH: f 1 = 21,395.85, f 1 = -1779.70 with nonzero SLCs. The v cr_lin stays the same when f 1 changes a lot in Fig. 7h. No GH point is observed within the calculated range.
It should be mentioned that, due to the symmetrical structure of the wheelset, the even order items are canceled out during force analysis. Thus, the effects of a 2 as well as c 2 , e 2 , and f 2 , described in Eqs. (1)-(4), on the bifurcation analysis are neglected. In case of uneven wear, which introduces differences on the even order items of the polynomials for the left and right wheel, the even order items should be included, while this case was not investigated in researches [4,9,16], and this work also.  a a 1 , b a 3 , c c 1 , d c 3 , e e 1 , f e 3 , g f 1 , and h f 3

Analysis of LPC bifurcation
For railway vehicles, one can take v cr_lin as an upper limit to determine the operating speed with a reasonable safety margin considered in the supercritical Hopf bifurcation cases. However, the hunting motion with a magnitude up to 10 mm may occur even when the running speed is less than v cr_lin in the subcritical Hopf bifurcation cases. Then the limit cycles are calculated through the continuation algorithm by Matcont [37] to determine the nonlinear critical speed (v cr_non ). LPC bifurcation curves and Hopf bifurcation curves are illustrated in Fig. 8, and a few points are listed in Table 3. LPC bifurcation appears only with a positive FLC, i.e., a subcritical bifurcation. To observe its relationship with each parameter, a 3 takes a value of 200 rather than the default value -499.40 except in Fig. 8b. Each figure can be divided into three parts by Hopf curve and LPC curve. The area above the Hopf curve is subjected to an unstable equilibrium surrounded by a stable limit circle with a magnitude close to 10 mm, which implies that the system would end up with a heavy hunting motion with a large magnitude. The area below the LPC curve is subjected to a stable equilibrium, which implies that the system would return to equilibrium. The area between them is subjected to a stable equilibrium surrounded by an unstable limit circle with a magnitude of 0-10 mm and a stable limit circle with a magnitude close to 10 mm, which implies that the system would return to the equilibrium if the disturbance is small, while it would end up with heavy hunting motion with a large magnitude if the disturbance is large. From Fig. 8a, one can tell that the v cr_non decreases with a ratio inversely proportional to a 1 as the v cr_lin . It drops by 40.34% when a 1 increases from 0.05 to 0.15.
Meanwhile, the v cr_lin and v cr_non get closer as a 1 increases. This shows coincidence with the decreasing FLC, which gets closer and closer to zero. When a 3 takes different values in Fig. 8b, the v cr_lin remains almost the same, but the v cr_non varied significantly as well as the FLC in Table 3. The larger a 3 , the larger the FLC, and the smaller the v cr_non .
In Fig. 8c, the v cr_non decreases significantly with c 1 ; thus, a maximum of c 1 should be imposed to ensure a reasonable v cr_non . The v cr_lin decreases little within the calculated range. Another GH point is observed at point: c 1 = 1140.8416. However, it is unlikely that c 1 can be so large in reality. As in Fig. 8d, the v cr_non decreases 0.1 m/s while c 3 changes from -1 9 10 5 to 1 9 10 5 .
As in Fig. 8e, g, the influence of e 1 and f 1 on v cr_non is alike. v cr_non rises with them sharply at first then linearly and finally reaches a GH point and coincides  with the Hopf curve. As in Fig. 8f, h, v cr_non changes a little while e 3 and f 3 vary a lot. Their influence on either the Hopf bifurcation or the LPC bifurcation can be ignored within the calculated range.

Stability analysis under wheel wear
Preceding studies show that geometry nonlinearities of wheel profile have a significant effect on the bifurcation diagrams of vehicle, in which the wheel profile is a standard one with smooth curvatures and radius, whereas this section concerns the effect of wheel wear on vehicle stability, which introduces highly unsmooth characteristics. The measured worn profiles on the field for a high-speed train operated in China Railway were employed for numerical investigations.
6.1 Geometry interpolation of worn profiles Figure 9a compares the measured wheel profiles between several running mileage (divided by 10,000 km), where the wear dominates at ± 15 mm around its origin of the coordinate system as depicted in Fig. 9b. In those cases, as in Fig. 10a, the rolling radius increases almost linearly with a ratio proportional to running mileage when the lateral displacement is larger than -2.5 mm. This is consistent with the fact that the equivalent conicity increases with running mileage, as depicted in Fig. 10b. From Fig. 10c, the contact angle varies the most with a worn profile of 0.184 million mileage. If a cubic polynomial is used to fit the function between contact angle and lateral displacement, the linear term could not vary that much as in Sects. 4 and 5.
As in Fig. 10d, the radius of curvature of a new wheel is much smoother than the others. For the worn profiles, it seems that sine curve fitting will be more proper.   In summary, the rolling radius of a new profile can be taken as a linear function of lateral displacement approximately, i.e., the equivalent conicity can describe the character. Its contact angle and curvature differ little, which has been ignored in dynamic simulation. However, the nonlinearity should be considered for worn profiles.

Hopf bifurcation analysis under wheel wear
The v cr_lin for each worn profile is listed in Table 4. It increases with running mileage at first then decreases. A negative FLC indicates a supercritical bifurcation while a positive one indicates subcritical bifurcation. Compared with the discussion in Sect. 5, it is unlikely that the terms of either the rolling radius and contact angle vary so much that the sign of FLC is changed. One can assume that the deviation of FLC is mainly caused by the radius of curvature of the wheel. However, its mathematical expression is not as accurate as wanted, a more reasonable fitting method would be studied in future work.
As in Fig. 11, the required suspension (steering) stiffness for a Hopf bifurcation of 70 m/s within the whole reprofile period varies with a maximum of 25%. Since the wheelset is connected to a fixed body (ground) through the suspension, the system would be more stable with larger steering stiffness. Larger stiffness would increase the running stability; however, it would worsen the wheel wear and safety on curved tracks. In addition, active control of suspension parameters can be introduced to solve the trade-off of suspension needs between the hunting stability and the curve negotiating behavior, while the bifurcation characteristics under the active control scheme have not been discussed yet.
However, the track irregularities have not been considered yet. Pieces of literature state that the excitation frequency and amplitude, as well as the suspension damping and its viscoelastic effect, have obvious influences on the system bifurcation performance and vibrating characteristics, see [38,39] for instance. It will be studied during future work.

Summary, conclusions and discussions
This investigation aims to study the stability and bifurcation characters of a railway wheelset under fully nonlinear wheel/rail interaction, including the contact geometry and forces. The geometry nonlinearities of wheel and rail profiles come from field measurements, including the rolling radius, contact angle, and curvatures, etc. To make the model more reliable, existing flange force models and their effects on the stability bifurcations are firstly compared. Then the geometry nonlinearities of profiles, namely the linear or cubic terms, on the Hopf and LPC bifurcations are thoroughly assessed. The following conclusions can be drawn, 1. The exponential flange model does not change either the Hopf bifurcation point or the bifurcation type and works well simulating the wheel flange in reality, while the others, polynomials, and reciprocal models, predict unreal bifurcations concerning reality. 2. The linear critical speed (Hopf bifurcation point) increases with the decreasing of a 1 (linear term of polynomial coefficients for fitting the rolling radius) significantly implies the equivalent conicity of wheel/rail contact, whereas c 1 (linear term of polynomial coefficients for fitting the contact angle) noticeably changes the Hopf bifurcation point, which decreases at first then increases with c 1 . But the linear critical speed varies slowly with e 1 and f 1 (linear term of polynomial coefficients for fitting the lateral radius of curvature of wheel and rail, respectively). 3. GH point observed with changing the value of a 1 , a 3 , c 1 , e 1 , and f 1 , i.e., they can change the bifurcation type; however, it is unsure that if then GH point exists in reality, while the cubic terms for the contact angles and the radius of curvature of wheel and rail profiles, c 3 , e 3 , and f 3 , do not change either the linear critical speed or the Hopf bifurcation type. 4. FLC can be used to predict the speed gap between the nonlinear critical speed and the linear one in the same model with one changing parameter. 5. The influence of e 1 and f 1 on the nonlinear critical speed is alike. It rises with them sharply at first then linearly and finally reaches a GH point and coincides with the Hopf curve. 6. Field measurements show that the equivalent conicity of wheel/rail contact increases approximately linearly with the increased running mileage, but jumps under certain cases. The shape of the conicity curve shall be concerned when performing bifurcation analysis rather than an equivalent constant commonly used in linearized and simplified models. The presented nonlinear model shows great priorities.
In summary, the cubic term for the contact angle and the radius of curvature of wheel and rail profiles can be ignored in the bifurcation analysis and dynamic calculations. However, it is crucial to take the linear term into consideration. The nonlinearity of rolling radius should be treated carefully since their influences on both Hopf bifurcation and LPC bifurcation are significant. An upper limit of a 1 , a 3 , c 1 and a lower limit of e 1 and f 1 should be imposed to ensure enough stability.
To solve the trade-off between the needs of suspension stiffness/damping for different degrees of wear of wheel, concerning the bifurcation characteristics under nonlinear wheel/rail contact, active primary/secondary suspensions are key techniques for high-speed trains to prolong the running mileage within a wheel profiling cycle, to save maintenance cost. While the bifurcation characteristics under the active control scheme have not been discussed yet, which will be studied in further work.