Dynamic Investigation of the Hunting Motion of a Railway Bogie in a Curved Track via Bifurcation Analysis

The main purpose of this paper is to analyze and compare the Hopf bifurcation behavior of a two-axle railway bogie and a dual wheelset in the presence of nonlinearities, which are yaw damping forces in the longitudinal suspension system and heuristic creep model of the wheel-rail contact including dead-zone clearance, while running on a curved track. Two-axle railway bogie and dual wheelset were modeled using 12-DOF and 8-DOF system with considering lateral, vertical, roll, and yaw motions. By utilizing Lyapunov’s indirect method, the critical hunting speeds related to these models are evaluated as track radius changes. Hunting defined as the lateral vibration of the wheelset with a large domain was characterized by a limit cycle-type oscillation behavior. Influence of the curved track radius on the lateral displacement of the leading wheelset was also investigated through 2D bifurcation diagram, which is employed in the design of a stablemodel. Frequency power spectra at critical speeds, which are related to the subcritical and supercritical bifurcations, were represented by comparing the two-axle bogie and dual wheelset model. The evaluated accuracy to predict the critical hunting speed is higher and the hunting frequency in unstable region is lower compared to the dual wheelset model.


Introduction
The prior issue, which should be clarified in railway vehicle design under rail load, is the dynamic response of the vehicle having clearances between wheel tread-rail, wheel flangerail contacting forces, and dry friction in suspension system. When the velocity of the railway vehicle increases, the vehicle evolves into a less stable state and expresses strictly oscillations due to the velocity dependency of "lateral and yaw coupled" stability, commonly called "hunting" [1]. The motion dynamics is highly nonlinear by the reason of the effects of nonconservative wheel-rail creep forces. To design safe and robust vehicle, it is essential to investigate the influence of the railway vehicle parameters on the hunting motion. Hunting is characterized by a limit cycle-type oscillation. The most important nonlinearities, which have influence parameters of hunting or chaotic motion occurrence, are classified such as clearances between the components, the wheel flange contact forces, and dry friction in suspension components [2]. The real cause of the hunting phenomenon is the tangential contact forces in combination with the longitudinal and lateral suspension forces. "Bifurcation" or "branching" is a phenomenon that takes place in nonlinear problems, when their solutions rely upon one or more parameters, which characterize the problem. As examples of such parameters, Reynolds number in fluid mechanics, the dimensionless load in buckling of structures, and the speed of the train in the problem of bogie oscillations can be given [3]. The stability assessment in railway industry necessitates the chance to apply the stability analysis on complex vehicle models under realistic operation conditions. Polach [4] has applied the use of nonlinear calculations and bifurcation analysis when dealing with the running stability during vehicle design and development in railway industry and discussed the relationship between the bifurcation diagram and the assessment of safety risk and the dynamic behavior. Wu and Chi [5] have formulated a suspended single wheelset model, which was adopted to the roller rig tests to investigate the effect of 2 Mathematical Problems in Engineering the vehicle parameters on the features of Hopf bifurcation. The lateral and yaw motions of the wheelset were considered in their model. Suspension elements and creep forces were taken into account as linear properties. To study the features of Hopf bifurcation of the wheelset model, the continuation method was used and examined the combined influence of parameters and coupling regions between the parameters. Kaas-Petersen [6] has analyzed the hunting motion in the Cooperrider bogie and has revealed that the behavior of the railway bogie has been strongly influenced by wheel flanges. If the bogie has wheel without having a flange, a symmetry breaking bifurcation is emerged with a transitive process, which is changing from a symmetric periodic oscillation state to an asymmetric periodic oscillation state. If the bogie has flanged wheels, which is represented as a dead-band spring, chaotic behavior appears [7]. Zboinski and Dusza [8] have investigated the nonlinear lateral stability of railway vehicles with different pairs of wheel/rail profiles in a curved track by using self-exciting vibration and bifurcation theories. Lieh and Haque [9] have studied the excited behavior of passenger and freight vehicles on a tangent track by the reason of harmonic variations in the conicity via linear models and the effect of primary and secondary stiffness on parametric excitation. Zhang and Dai [10] have built a lateral mathematical model of railway wheelset to study the effect of yaw damper on the stability under the rigid constraint between the wheelset and bogie. Center manifold theorem, the method of normal form, and Poincare method have been used to reduce the model to a planar dynamical system and which kind of bifurcation will occur at the critical speed has been determined. Ahmadian and Yang [11] have searched on the analytical representation of Hopf bifurcation and hunting behavior of a rail wheelset with nonlinear primary yaw dampers and wheel-rail contact forces. They have proved that the nonlinearities in the primary suspension and flange contact contribute significantly to the hunting phenomenon. The critical speed and the nature of bifurcation are influenced by the nonlinear elements. Reza [12] has investigated the analytical formulation of bifurcation, nonlinear lateral stability, and hunting behavior of rail vehicles in a tangent track by using Bogoliubov method. This model includes also a nonlinear primary yaw dampers and flange contact. Due to the three-dimensional creep forces acting between the wheel and rail, the lateral movement of the wheelset can be oscillated around the equilibrium position with increasing or decreasing amplitude over time, even in the case of positive conicity. Factors such as wheelset yaw and track irregularity also reinforce this kind of oscillation. The resulting vibration must remain acceptable to provide certain comfort and safety requirements. Depending on the design of the railway vehicle, it may experience severe oscillations, which may be the cause of the source of disturbance or derailment above a certain operating speed. The speed at which the railway vehicle becomes unstable is called "critical speed." Above the critical speed, the vehicle is subjected to much higher forces due to the increasing oscillating motion and to the collision between wheel flange and the rail. The impact between the flange and the rail is known as the secondary contact point. The contact point between the wheel tread and rail head is also referred to as the primary contact point [13]. Two-point contact occurs in the outer wheel of the leading wheelset of a rail vehicle moving on a curve: one contact patch among the wheel tread and the rail head and the other between the flange and the rail gauge. The two-point contact situation is necessary in the analysis of hunting and derailment. Piotrowski and Chollet [14] described the force distribution over the tread and flange contact zones for a rolling wheelset. The contact forces for both the contact zones are determined using Kalker's empirical nonlinear theory under the condition of a rigid wheelset [15]. Linear and nonlinear critical speeds were obtained while concentrating the change of the parameter effect in differing the speeds with altered lateral stiffness in primary suspension system. Scheffel [16] analyzed the combined effect of the longitudinal and lateral suspensions on the hunting stability in a four-wheeler model. According to this study, the bogie hunting can be restrained if the speed, where kinematic frequency of the wheelset accompanies the natural frequency of the bogie vibration, is low compared to the critical speed of wheelset hunting. True et al. [17] found that the allowed speed on the curved track and the critical speed on a tangent track can be higher than the critical speed only in curves by investigating the quasi-stationary dynamics of the Cooperrider bogie running on a curved track with a realistic wheel-rail contact geometric relation. The excitation terms in the equations of motion that was derived for stability analysis were validated by using a roller rig model. A comprehensive approach to the dynamic response to large discrete inputs, including both suspension and wheel-rail contact nonlinearities, was conducted with a full-scale experimental validation. It can be stated that irregularities are perturbed continuously along the track. The approach proposed by stochastic process theory was applied to the lateral motions of a restrained wheelset. Full railway vehicle responses were established as an indication of ride quality and for passenger comfort assessments. The first measurements of power spectral density using a specially developed trolley-based measurement system were carried out. Extensive measurements of the power spectra of irregularities of the track have been made on railways resulting in inputs utilized in design [18]. There is a phase difference between the inputs of the front and rear wheelsets because the rear wheelsets pass the same portion of track after the front wheelsets. If the same input to the rear wheelsets with a time delay was given, the phase difference can be modeled. The rail joints on the right and left tracks may be staggered. This is another source of a phase difference, which causes roll motion of the vehicle. The effects of staggered rail joints can not be captured by using the 10-DOF model in which the roll degree of freedom is not included [19]. The nonlinear dynamical characteristics of a railway vehicle are checked into thoroughly by applying two different wheelrail contact models: a heuristic nonlinear friction creepage model derived by using Kalker's theory and Polach model including dead-zone clearance. Two models are matched with the quasi-static form of the LuGre model to obtain a more realistic wheel-rail contact model. LuGre model parameters are determined using nonlinear optimization method, whose objective is to minimize the error between the output of the Polach and Kalker model and the quasi-static LuGre model for specific operation conditions, in our previous research [20]. However, in the case of high values of spin, the heuristic theory leads to unsatisfactory results [21].
In this paper, a low frequency bogie and two-axle wheelset's hunting problem of a conventional high speed railway vehicle were examined with respect to the various radius of the curved track. The models have been constructed by considering all of the possible design concerns from the viewpoint of hunting stability which existed in the literature to estimate the critical hunting speed more accurately. A complex nonlinear wheel-rail contact relation is applied in our analysis of the bifurcation for the lateral displacement of the leading wheelset with respect to the different radius of the curved track. The results evaluated from the heuristic model are in good agreement with Kalker's simplified theory. The wheel flange force is modeled as a linear spring with a nonlinear damping including a dead zone caused by the wheel-rail clearance. This investigation emphasizes the influences of the system suspension nonlinearities and the wheel-rail interface nonlinearities on Hopf bifurcation. Hunting frequencies and critical velocities are found by using PSD (power spectral density) and Lyapunov's indirect method for bogie and dual wheelset model. The relationship between the hunting speed and radius of the curved track parameters of the system is introduced. Two-dimensional bifurcation diagrams to study the Hopf bifurcation in the system are demonstrated. It was deduced that the railway bogie is getting destabilized when the kinematic frequency of wheelset approaches to a natural frequency of the bogie on the suspension system. From the knowledge acquired from this investigation, one concludes that the critical hunting speed information evaluated from two-axle bogie model is more accurate as compared with the dual wheelset model. The acquired results are suitable for railway industry from the point of view of experimental and operational systems.
The remainder of this paper is arranged as follows. Section 2 contains the formulations of the equations of motion of a two-axle bogie and dual wheelset, respectively. Numerical analysis of the hunting behavior of these systems is done in Section 3. Results and discussion are presented in Section 4 and concluding remarks are given in the final section. Figure 1 is a schematic of the dual wheelset and bogie system, where stiffness and damping elements have been incorporated into the suspensions. Figure 2 represents the free-body diagram of the wheelset. Dual wheelset model is considered to be 8-DOF system and bogie model is considered to be 12-DOF system. The difference between the two-axle bogie and dual wheelset model is that the lateral, vertical, rolling, and yawing motions of the bogie frame are not considered in the dual wheelset model; therefore the lateral, vertical, yaw, and roll motions of the bogie frame are not coupled to the dual wheelset model. It means that their coordinates are equal to zero when suspension forces are calculated. The bogie frames and wheelsets are connected through the suspensions.

Description of the System.
Suspensions are modeled with three spring-damper systems that are capable of movement in the longitudinal, lateral, and vertical directions. Lateral, vertical, yaw, and roll motions of the system components (bogie frame and wheelsets) are coupled through these suspensions. The contact mechanisms between the wheel and rail surfaces are complex and are of great importance because they partially govern the dynamic behavior and hunting stability. The contacts between the wheel and rail can be classified into primary and secondary contacts that occur at the interfacial surface and at the wheel flange and rail side, respectively. The primary contacts induce nonconservative creep forces in the elliptical contact region. The secondary contacts impact the rail side when the actual clearance between the wheel flange and the rail is zero. In this study, the wheel is assumed to be conical, and the rail is modeled as possessing a knife-edge and having known lateral and vertical stiffness. This kind of linear geometry can be compensated with the stiffness force representing rail/wheel flange contact. Such action is not adequate for classical twopoint contact (wheel tread/rail, wheel flange/rail). Such type of contact is more often expected in tramway systems where curve radii are much smaller, the profiles can be different than railway ones, wheel and rail wear are usually higher, and maintenance is usually weaker than in railway systems. For real nonlinear wheel and rail profiles geometry, influence of the flange area starts gradually and relatively mildly with the fluent increase of contact angle corresponding to bigger lateral wheel/rail shifts. Such a simplification is preferred for ease of calculation (computational effort). It is claimed that it is compensable because it compares two identical systems (two-axle bogie and dual wheelset).

Wheel-Rail Contact Model
. A heuristic contact model, in which the nonlinear effect of the adhesion limit is incorporated, is derived to reveal the creep mechanism. Kalker's linear creep theory has been modified with the function of creep force saturation to evaluate the creep forces acting on the wheel-rail interface [24,25]. The unlimited creep forces and moments at the left and right contact planes are solved by utilizing Kalker's linear theory as * = − 33 K px , C px Figure 1: Schematic of the two-axle railway bogie.
where the subscripts " " ( = 1 for the leading wheelset and = 2 for the trailing wheelset) show the positions of the wheelsets in the railway vehicle and the subscripts and represent the left and right hand sides, respectively. The saturation constant is calculated by using the unlimited resultant creep force as where The vertical and lateral components of normal forces, which affect the primary contact patch, are determined with the help of the wheelset configuration information as follows: The vertical reaction force acting on each wheel is modeled as a spring force, which is an implicit function of the axle load ( axle = ext + ), suspension, and inertial forces of the wheelset. The flange contact force at the secondary contact patch is evaluated as a spring reaction force as The logistic functions used in (2) and (5) are utilized to characterize the nonlinear dead-band behaviors of the vehicle originated from the saturation of the creep constant and the flange contact without falling any numerical singularity. With the help of the numerical analysis and the actual contact mechanism, the saturation constant and flange contact force proposed are converged to their real values by fusing hyperbolic tangents into the functions.

Derivation of the Equations of Motion.
Three sets of Cartesian coordinate system are depicted in Figure 3. The coordinate system " , , " has its origin at the track center line and moves at a constant forward velocity with respect to a fixed inertial reference frame. The coordinate system " , , " is an intermediate frame that is rotated through an angle about the axis . The axes " , , " make up the wheel axle set body coordinate system, which has its origin at the center of mass [24].
The equations of motion for the bogie frame and wheelsets in the lateral, yaw, vertical, and roll directions are based on the following equations [22], respectively.̈= The moments generated by the nonlinear longitudinal damping forces ( ) are included in vertical suspension moment ( ). These damping forces can be represented by fourthorder polynomials of [11].
where =̇is the relative longitudinal velocity between the bogies and their wheelsets (which leads to the yaw motion). The additive moments generated by these forces are represented as = The nonlinear heuristic creep forces and moments can be represented by their corresponding linear creep force multiplied by the associated saturation factor such that The position vectors can be formed by using geometric configuration of the vehicle-rail system as The suspension forces and moments acting on wheelsets and bogie, which are induced from the lateral stiffness and yaw damping of the primary suspension, can be expressed as follows:

Lyapunov's Indirect Method to Determine the Critical
Speed. This kind of equations of motion has nonlinear nature. In this method, to calculate the critical hunting speed, the dynamical equations of motion should be expressed in the state space form aṡ( where ( ) denotes 16-dimensional state vector for wheelset case and 24-dimensional state vector for bogie case. [ ( )] is also same dimensional nonlinear vector, that is, a function of the state variables [27,28]. To obtain the equilibrium points, the nonlinear algebraic equations (6) must be solved. Taylor series expansion is applied to the right side of (13) at the equilibrium point. The linearized version can be obtained by utilizing the following matrix/vector form aṡ where ≡ [ ( )] = and is the Jacobian matrix calculated at the equilibrium point from the linearized system matrix.
With the help of the evaluated eigenvalues of the system, the hunting stability problem can be clarified.
It can be inferred that when all eigenvalues of matrix " " are placed in the left-half complex plane, the equilibrium point, which corresponds to the actual nonlinear system, is asymptotically stable. When all eigenvalues of matrix " " are placed in the left-half complex plane except at least one of them being on the imaginary axis, one can not make a judgement of stability from the linear approximation due to the marginally stable case. When at least one eigenvalue of matrix " " is strictly in the right-half complex plane, the equilibrium point is unstable for the nonlinear system [29].

Definition of the Limit Cycle-Type Oscillation.
Oscillation is one of the most important phenomena that occur in dynamical systems. A system vibrates when it has a nontrivial periodic solution.
for some > 0. "Nontrivial" term is used to exclude constant solutions corresponding to equilibrium points. The image of a periodic solution in the phase portrait is a closed trajectory, which is usually called a periodic orbit or a closed orbit [30]. The oscillation of nonlinear systems due to a center equilibrium point is not structurally stable. Infinitesimally perturbations may change the type of the equilibrium point to a stable focus (decaying oscillation) or unstable focus (growing oscillation).

Hopf Bifurcation Characteristics of Two-Axle Railway
Bogie. Bifurcation is defined as a change in the equilibrium points or periodic orbits or in their stability properties, when a parameter varies. Hopf bifurcation points out a critical parameter value, by which the equilibrium point of the system becomes unstable and produces the limit cycle behavior.
The Hopf bifurcation occurs generally in the supercritical and subcritical forms as shown in Figure 4.
Supercritical bifurcation indicates the limit cycle and the stable equilibrium point of the system divided into two sides of the Hopf bifurcation. Controversially, the subcritical bifurcation of the system refers to the limit cycle and a stable equilibrium point found on the same side of the Hopf bifurcation. The solid lines show a stable limit cycle, and the dotted lines show a nonstable limit cycle. The speed at point , which is named Hopf bifurcation point, is the linear critical velocity. is the saddle-node bifurcation point. The speed at point is the nonlinear critical speed. The limit cycle-type oscillations are generated when the vehicle speed exceeds the linear critical velocity. If the critical speed at this point is exceeded, the maximum value of the limit cycle amplitude increases or decreases progressively as speed increases or decreases [31,32]. In the subcritical bifurcation case, while vehicle speed is left of the point , the system stands absolutely in stable state. The system exhibits hysteresis phenomenon. When the speed exceeds point , the system suddenly produces significantly limit cycle and jumps to  large stable limit cycle or periodic oscillation. If the speed reduces until point , the periodic oscillation damped to stable equilibrium position [33]. If the value of the speed is between and point, one cannot decide about the characteristic of the system stability. The system tends to instability when the perturbation is large, by contrast with the perturbation which is small; then the system remains in stable condition. Stability can be regained only when the vehicle speed is reduced to [34].

Methods of Computation.
Several methods exist for solving nonlinear multibodied equations. Two common methods are Newton Raphson's method and the Quasi-Newton method. Newton Raphson method is a numerical method for solving simultaneous nonlinear equations. Quadratic convergence is supplied by this method. The algorithm for implementing this method is where −1 is inverse Jacobian matrix of ( ), is initial conditions for iterations, and +1 is the next state. Newton Raphson algorithm terminates only when the function ( ) is close to zero. Newton Raphson's method is less efficient than Quasi-Newton method. In Newton Raphson's method, the Jacobian matrix has to be evaluated for every iteration but with the Quasi-Newton method, a single Jacobian matrix is determined and thus used in the iteration. The algorithms can be improved by implementing Gauss-Newton, Levenberg-Marquardt, Trust-Region-Reflective method, and so forth. For example, Levenberg-Marquardt algorithm is a numerical optimization method used to reveal the solution of nonlinear function by minimization over a region of parameters of that function [35].
Runge-Kutta 4th-order method is a numerical technique used to solve ordinary differential equation of the form = ( , ) , This method is utilized to solve high order ordinary differential equations or coupled (simultaneous) differential equations. In the forward Euler method, the information on the slope or the derivative of at the given time step to extrapolate the solution to the next time step was used. Runge-Kutta methods are a class of methods which judiciously utilizes the information on the "slope" at more than one point to extrapolate the solution to the future time step. In contrast to multistep methods, the Runge-Kutta method, Mathematical Problems in Engineering

Determination of Critical Hunting Speed with respect to Various Radius of the Curved Track.
Lyapunov's indirect method was performed on the linearized system to investigate the variation of the critical hunting speed with respect to the different radius of the curved track parameters, which are changed between 250 and 10000 m. The equilibrium points were obtained by solving system of nonlinear equation with Levenberg-Marquardt algorithm iteratively. In order to solve (13), with zero initial condition, Runge-Kutta method of order four and the parameter values in " Table 1" were used. The 8-DOF dual wheelset and 12-DOF bogie model that considers flange contact, the heuristic creep, and nonlinear damping, were employed. The simulations were performed on Matlab&Simulink. The results are given in Figure 5 for dual wheelset and bogie model.
Comparing critical hunting speeds obtained from Lyapunov's indirect method and from numerical simulation, the relative error is approximately %8,75. The stability of the original nonlinear system then cannot be determined exactly using the linearized model.

Analysis of the Limit Cycle Motion.
At the hunting speed, a self-sustained vibration of the vehicle occurs, and the resulting motions show the limit cycles in the phase planes. Numerical simulations are performed with the linear elastic rail model ( = 1) and the radius of the curved track value = 6250 m. The lateral displacement and yaw rotation phase portraits and the time dependence of lateral displacement for , above and below the critical value cr-NL , are illustrated in Figures 6-9 both in bogie and in dual wheelset model, respectively.
One can conclude from below the critical hunting speed that the vibrations tend to the equilibrium point and stable motion reaches for freely long time behavior. By increasing velocity, the vibrations amplitude increases. Above the critical hunting speed, the vibrations of the limit cycle motion are preserved and the vibration amplitudes increase for freely long time.

Bifurcations of the Leading Wheelset.
The wheelset and bogie considered in this study include several nonlinear elements, such as the heuristic creep, flange contact, and nonlinear damping. The dynamic behavior of such a nonlinear system is highly dependent upon its initial state. As the forward velocity increases, the vehicle becomes less stable and ultimately exhibits rigorous oscillations [1]. In the bifurcation theory, appearance of hunting is represented by Hopf bifurcation of a fixed point. The motion is stable below a certain "critical" forward velocity (or exhibits a stable fixed point in phase space). Above the critical speed, the hunting appears as an undamped or damped vehicle motion through the tangential forces that depend on velocities through the creepages relations (the fixed point loses its stability and a limit cycle bifurcates from it). The decay rate of the phase portrait of wheelset depends on control parameter. Supercritical Hopf bifurcation occurs when a stable spiral changes into an unstable spiral surrounded by a small limit cycle [39]. The stability and dynamic behavior of the nonlinear system can exhibit divergent result depending on its initial conditions [40]. In order to analyze the influence of the radius of the curved track parameter on the hunting behavior, the forward speed V varies between 10 and 180 m/s. The stable response of lateral displacements versus speed with respect to different radius of the curved track values is plotted for 5 sec. Bifurcation diagrams are drawn for bogie and dual wheelset model at range of curved track radius of 250 m, 500 m, 750 m, 1000 m, 2500 m, 6000 m, and 10000 m, respectively, in Figures 10 and 11.
It can be derived from the figures that the maximum lateral displacement is getting larger when the radius of the curved track decreases and the critical hunting speeds occurred at lower speed while the radius of the curved track   [1,22,23].

Parameters
Numerical value Vehicle body, bogie frame, and wheelset masses (kg) V = 34000, = 3000, = 1800 Nonlinear longitudinal damping coefficients (kNs/m) 1  decreases. It should be noted that use of the fully nonlinear models could lead to results different than those shown in both figures. For example, subcritical properties might eventually be reviled.

Power Spectral Density Approach.
To specify power spectrum of the lateral displacement of the leading wheelset, power spectral density approach should be carried out [23]. The result of this analysis is presented for subcritical hunting speed in Figure 12 and for supercritical hunting speed in Figure 13. As seen in Figure 12  ( cr-NL < < cr-L ) term referred to in Figure 12 states the speed between the linear and nonlinear critical speed. Frequency response for the supercritical hunting speed is illustrated in Figure 13. In Figure 13, main frequencies for the bogie and dual wheelset model are 2.734 and 2.78 Hz, respectively. RMS values for bogie and dual wheelset model are 77.9542 and 74.5568, respectively. The hunting frequency calculated via the dual wheelset model is greater than that derived using the bogie model. The railway bogie loses its stability when the kinematic frequency of wheelset approaches to a natural frequency of the bogie on the suspension system.
The individual wheelset kinematic mode can emerge if its frequency is higher than the bogie frame natural frequencies [41]. The frequency of the bogie mode is the kinematic frequency of the coupled wheelset mode. At low speeds, this kinematic frequency is lower than the bogies natural frequencies, which are usually around 1 Hz. The motion of the vehicle is primarily determined by the wheel/rail forces and by the primary suspension forces. The car body moves without affecting the motion of the vehicle significantly. If the speed of the vehicle is gradually increased, a special situation occurs in the kinematic frequency of the vehicle. Therefore, the car body hunting problem can usually be solved by increasing the secondary lateral damping ratio. A well designed vehicle can experience increased but stable motions near the speed at which the wheelset kinematic frequency equals the bogie's natural frequencies.

Conclusion
The system dynamic of the two-axle bogie and dual wheelset was modeled using 12-DOF and 8-DOF system together with the lateral, vertical, yaw, and roll displacements, respectively. Phase portraits were drawn to analyze the railway bogie movement and related limit cycle motion. Runge-Kutta order of four methods was used to solve the nonlinear equations of motion. Hopf bifurcation in two-axle bogie and dual wheelset was also studied through a nonlinear model, which includes nonlinear wheel-rail contact and nonlinear yaw damper suspension. Critical hunting speeds were evaluated with respect to the different radius of the curved track through Lyapunov's indirect method and the results were compared with the simulations in Matlab&Simulink. Long term behavior exhibits that the oscillations tend to the equilibrium state for speeds lower than the critical hunting speed. Stable movement is ensured lower than the critical speed. When the railway bogie is at critical hunting speed, the limit cycle motion is raised. Above the critical speed, the vibrations amplitude increases. The hunting frequencies of the lateral displacement of the wheelset were calculated via PSD method for two different models. In consequence of the simulations, it can be concluded that when the curved  track radius decreases, the critical speed reduces. The critical hunting speed and hunting frequencies obtained from dual wheelset model are higher than those derived from twoaxle bogie model. It can be deduced from the analysis done that the critical speed obtained from two-axle bogie model becomes more accurate.    Left and right wheel of the rolling radii, respectively , , , : Lateral, vertical displacements and yaw, and roll angular displacements of the wheelsets, respectivelẏ ,̇,̇,̇: Lateral, vertical velocities and yaw, and roll angular velocities of the wheelsets, respectively * , * :

Nomenclature
Linear creep forces of the left wheels in the longitudinal and lateral directions given by Kalker's linear theory, respectively * , * : Linear creep forces of the right wheels in the longitudinal and lateral directions given by Kalker's linear theory, respectively , : Contact angle of the left and right wheels, respectively * : Linear creep moment of the left wheels in the vertical direction given by Kalker's linear theory * : Linear creep moment of the right wheels in the vertical direction given by Kalker's linear theory : Saturation constant in the heuristic creep model : Nonlinearity in the heuristic creep model , : Nonlinearities of the left and right wheels in the heuristic creep model, respectively , : Components of normal forces on the left wheels in the lateral and vertical directions, respectively , : Components of normal forces on the right wheels in the lateral and vertical directions, respectively , : Resultant normal forces on the left and right wheels on the equilibrium axis, respectively : Coefficient of friction : Wheel conicity : Flange clearance , : Lateral and vertical rail stiffness, respectively : Gravitational constant : C a n t a n g l e : W h e e l s e t m a s s Linear creep forces at the equilibrium coordinate system of the left wheels in the longitudinal, lateral, and vertical directions, respectively , , : Linear creep forces at the equilibrium coordinate system of the right wheels in the longitudinal, lateral, and vertical directions, respectively , , : Linear creep moments at the equilibrium coordinate system of the left wheels in the longitudinal, lateral, and vertical directions, respectively , , : Linear creep moments at the equilibrium coordinate system of the right wheels in the longitudinal, lateral, and vertical directions, respectively , : Nonlinear heuristic creep forces of the left wheels in the longitudinal and lateral directions, respectively , : Nonlinear heuristic creep forces of the right wheels in the longitudinal and lateral directions, respectively , : Nonlinear heuristic creep moment of the left and right wheels in the vertical direction, respectively , , : Longitudinal, lateral, and vertical components of the contact position vector on the right wheels, respectively , , : Longitudinal, lateral, and vertical components of the contact position vector on the left wheels, respectively , : Lateral and vertical suspension forces of the wheelsets, respectively , , : Suspension moments of the wheelsets in the longitudinal, lateral, and vertical directions, respectively , , : Roll, pitch, and yaw moment of inertia of the wheelsets, respectively , : Roll and yaw moments of inertia of the bogie, respectively , , : Primary longitudinal, lateral, and vertical stiffness, respectively , , : Primary longitudinal, lateral, and vertical damping coefficients, respectively 1 , 3 : H a l f o f t h e p r i m a r y l o n g i t u d i n a l a n d vertical spring arms 2 , 4 : H a l f o f t h e p r i m a r y l o n g i t u d i n a l a n d vertical damper arms ℎ : Height of the bogie mass center above the wheelset mass center 1 , 2 : Half of the primary lateral spring and damper arms : Distance between the wheelset and the bogie frame mass center.