In this section, the dynamic and quasi-static models are presented for the pummeling analysis approach. More specifically, the section is organized as follows: the dynamics models of the wagon using the commercial software SIMPACK® and VAMPIRE® are given in Sec. (2.1). In Sec. (2.2) we discuss the quasi-static model employed by the software GIRAFFE. Finally, in Sec. (2.3) the pummeling analysis for the case of a heavy-haul railway in Brazil is presented.
Dynamic simulations
The heavy haul wagon GDE-Ride Control was modeled in SIMPACK® and VAMPIRE®. It is a meter gauge three-piece bogie wagon with 27.5 tons per axel, employed for ore transportation in a Brazilian railway. The model of the wagon was validated by comparison with previous results from NUCARS® and experiments performed under supervision of the ore transportation company. Figure 2 shows the comparison of the results from the three software for a specific sector of the railway.
The primary suspension includes a shear pad between the wheelset bearing and the side frame pedestal. The pad adapter is made of an elastomeric material positioned between two steel plates [9]. It was modeled in SIMPACK® using three springs between the wheel axle and the side frame. The inner and outer springs feature a vertical stiffness and a vertical damping parameter (𝑧). The central spring has stiffness and damping in the longitudinal (𝑥), lateral (𝑦), and rotational (yaw -about 𝑧) directions. The magnitudes of stiffness and damping adopted include the influence of the pad. In VAMPIRE®, the primary suspension is defined by bushing elements, that have linear stiffness and damping, and non-linear stiffness elements. These elements together represent pitch and yaw stiffness and vertical, roll and yaw damping.
The secondary suspension is modeled in SIMPACK® with bushing springs. The spring system is composed of four springs arranged in a square distribution, with a vertical stiffness and a vertical damping parameter (z). A fifth spring is in the center, which has stiffness properties and damping parameters in all directions. In the same suspension system, the wedges with constant damping were modeled using discrete contact elements and considering dry friction between the contact surfaces. In VAMPIRE®, the secondary suspension is modeled as a seven-spring package with vertical, lateral and longitudinal linear stiffness and a central element of non-linear vertical stiffness. The wedge is represented by two elements of lateral non-linear stiffness, four friction elements and a spring with longitudinal linear stiffness.
The connection points between the wagon box and the bogies are modeled on the center plate and the roller side bearing. In SIMPACK®, dry friction elements were adopted on the side of the center plate for tangential contact and universal point contact elements for the vertical direction. For the contact model of the roller side bearing, a spring (bushing type) was used with the vertical stiffness of contact with a gap, since the type of balance support modeled is not of constant contact, like the one adopted for modeling the primary suspension. In VAMPIRE® the center plate is modeled by four damping elements and twelve friction elements at the ends of the plate, in addition to four non-linear stiffness elements at the same points. In the plate’s center, there are two other non-linear rigidity elements, longitudinal and lateral. The roller side bearing is represented by a non-linear vertical stiffness element with a gap of 8 mm.
As already mentioned, the loaded axle weight for this GDE wagon is 27.5 tonnes. Because of that weight, the wagon is considered a heavy-haul freight car. The models of the bogies for this wagon on SIMPACK® and VAMPIRE® are shown in Fig. 3 and Fig. 4, respectively.
The wheel-rail contact forces in SIMPACK® and VAMPIRE® are calculated from a simplified Kalker’s theory, which defines the relationship between the creepages and the creep forces. The Kalker’s linear theory (Eq. 1) uses the well-known Kalker coefficients (Cij) to determine the relationship between creepages and creep forces and moments. In this case, the longitudinal and lateral forces are given by a linear function [10],
$$\left[\begin{array}{c}{F}_{x}\\ {F}_{y}\\ {M}_{z}\end{array}\right]=-Gab\left[\begin{array}{ccc}{C}_{11}& 0& 0\\ 0& {C}_{22}& \sqrt{ab}{C}_{23}\\ 0& -\sqrt{ab}{C}_{23}& ab{C}_{33}\end{array}\right]\left[\begin{array}{c}{\xi }_{x}\\ {\xi }_{y}\\ \phi \end{array}\right]$$
1
.
F x , Fy and Mz are the longitudinal creep force, lateral creep force and spin creep moment. Cij (i,j = 1,2 and 3) are constants that depend on the Poisson coefficient ν and the ratio of the semi-axes of the contact ellipse patch. G is the modulus of rigidity, while a and b are the major and minor semi-axis of the contact ellipse, respectively. ξx, ξx and φ are the longitudinal creep, lateral creepages, and spin creepage. The full Kalker’s theory is implemented in the CONTACT software (software that numerically solves normal and tangential contact problems).
Normal and tangential W/R contact forces for the contact points are calculated using FASTSIM in SIMPACK®. It is a very fast method for calculating creep forces, widely used in contact W/R programs [11]. FASTSIM is a simplified approach based on that proposed by Kalker [12], in which the surface displacements are linearly related to the surface traction (Eq. 2) [13],
,
where p is the local tension and L is the flexibility constant, which is derived from the linear theory constants (Eq. 1).
For FASTSIM, Kalker proposed that L = Lx, Ly e Lφ. Thus, the flexibility constant is given by [14]:
$$L=\frac{\left(\left|{\xi }_{x}\right|{L}_{x}+\left|{\xi }_{y}\right|{L}_{y}+c\left|{\xi }_{\phi }\right|{L}_{\phi }\right)}{\left({\xi }_{x}^{2}+{\xi }_{y}^{2}+{c\xi }_{\phi }^{2}\right)}$$
3
$${L}_{x}=\frac{8a}{3{C}_{11}G}, {L}_{y}=\frac{8a}{3{C}_{22}G}, {L}_{\phi }=\frac{\pi a\sqrt{a/b}}{4{C}_{23}G},c=\sqrt{ab}$$
4
In SIMPACK®, the FASTSIM algorithm assumes an elliptical contact patch and discretizes it into a grid of mx x my elements (Fig. 5) [15]. Tangential forces are determined by integrating tangential stresses in grid elements numerically. Tangential stresses accumulate from the contact patch's front boundary until the local saturation limit is reached. Grid discretization in longitudinal direction is constant. The grid is automatically refined in the lateral direction around the spin pole if the pole lies within the contact patch. According to the contact patch width, the boundary indicates the minimum allowed grid width. It is determined from the initial discretization by my = mx: t = 1/[5(my − 1)]. SIMPACK suggest t = 0.02 for my = 11 as default, which is Kalker's recommendation [12, 16].
VAMPIRE® uses Hertzian contact theory to determine the normal forces and an off-line look-up table (LUT) to calculate the tangential contact forces. The use of LUT in the W/R contact requires a pre-calculated table that stores the values of the function determined with the CONTACT software, such as the normal and creep forces, in a set of points that constitute the domain of the table. Then, during a simulation, for a given input such as the relative position and orientation between rail and wheel, the output parameters required for the dynamic analysis are computed through the interpolation of data from the LUT. Thus, the online calculation of a cumbersome function, which is generally slow but accurate, can be promptly approximated through the linear interpolation of tabulated data [17].
VAMPIRE creates a multidimensional look-up table from the cross-sectional profiles of the wheels and rails using the Contact Data Generation preprocessor [18]. The program calculates the lateral displacement of each wheel relative to the rail at each time step. Based on this value, the contact data for each patch between wheel and rail will be looked up. The creepages for each wheel-rail contact patch are calculated using the position, rolling radius, and contact angle obtained for each contact patch as follows.
| \({\xi }_{x}=1\mp \frac{{l}_{0}.\dot{W}}{V}+\frac{r.\dot{P}}{V}+\frac{{l}_{0}}{R}+\frac{\dot{X}}{V}\) | (5) |
| \({\xi }_{y}=\left(\frac{r.\dot{P.W}}{V}+\frac{\dot{Y}}{V}\right)\text{s}\text{e}\text{c}\left(\delta \right)\) | (6) |
| \(\phi =\mp \left(\frac{\dot{P}}{V}\right)\text{s}\text{i}\text{n}\left(\delta \right)+\left(\frac{\dot{W}}{V}\right)\text{c}\text{o}\text{s}\left(\delta \right)\) | (7) |
In equations (5), (6) and (7), V is the vehicle speed; X,Y,W,P are longitudinal, lateral, yaw and pitch displacements of the wheelset in track axes, respectively; r are rolling radii of wheels; R is the curve radius of the track; l are the distances from the center of the wheelset to the left and right contact patch and δ is the contact angle.
Calculated creepages are then converted to non-dimensional values and used to interpolate non-dimensional creep forces from the Kalker look-up table. Vampire's default Kalker table is based on dry, uncontaminated wheel-rail contact conditions. In the run file, one can modify the coefficient of friction and Kalker factors to approximate other wheel rail contact conditions. As a result, non-dimensional creep forces are converted to actual creep forces at each point of contact between the rail and the wheel [18].
Quasi-static model
Higa et al. [8] proposed a simple 2D model to obtain the wheelset’s steady-state lateral position on a given track scenario, instead of solving a multibody transient dynamic simulation of the complete wagon traveling on a track. The model evaluates the normal contact interaction of a single wheelset with a single track section, with a fixed angle of attack for the wheelset. For each wheel/rail pair, the model can deal with multiple pointwise contacts. Its simplification premise is that the wheelset’s lateral steady state is governed by the quasi-static equilibrium between gravitational and centrifugal forces (active forces) and normal contact forces on the wheel/rail pairs (reactive forces), as shown in Fig. 6.
In this model the track is fixed, and the wheelset is free to move on the track’s transversal plane. Hence, the wheelset has three degrees of freedom, \({u}_{x}\) (lateral), \({u}_{y}\) (vertical) and \(\phi\) (rotation around longitudinal axis), and its position is given by the vector \(\varvec{d}\).
$$\varvec{d}=\left[\begin{array}{c}{u}_{x}\\ {u}_{y}\\ \phi \end{array}\right]$$
8
The track is composed of two rails positioned according to the rail cant angle, track gauge and superelevation. The wheelset is comprised of two wheels rigidly linked together and spaced by the back-to-back gauge. Both rails and wheels are modeled as rigid contact surfaces constructed respectively by extrusion and revolution of planar curves, meaning from the rails and wheels profiles.
For a tangent track scenario, only the gravitational load representing the entire wagon weight is applied. For a curved track, it is also applied the centrifugal inertial force dependent of the wagon mass, travel speed and curve radius. Since the model considers a single wheelset, the loads applied are divided proportionally to the number of wheelsets, four for the freight wagon studied here. To capture any effect of lateral load transfer between wheels, the loads are applied not at the center of the wheelset but at a node rigidly attached to the wheelset at a height equivalent to the entire wagon center of mass.
The model solution comprises the wheelset position, which gives the normal contact forces on each wheel/rail pair that will balance the applied loads. To obtain these contact points and normal forces, the model uses the master-master contact formulation [19], a penetration-based (soft contact) approach that can handle multiple contact points. Since the problem of finding multiple contact points and its normal forces is highly nonlinear, the model is solved adopting a dynamic relaxation procedure in which the nonlinear equations are integrated over time using Newmark implicit numerical integration to find the desirable equilibrium solution. The wheelset starts at a position above the rails; the loads are gradually applied; the wheels contact the rails giving a first solution for the normal contact forces; then the wheelset moves laterally seeking its equilibrium position automatically.
On the general 3D master-master contact formulation [19], both surfaces in imminent contact are treated as continuous surfaces with no discretization to pre-defined points; that means there is not a slave surface, as usual in node-surface formulations. The surfaces (\({{\Gamma }}_{A}\) e \({{\Gamma }}_{B}\)) are parametrized by their own convective coordinates (\(\zeta ,\theta\)) and the contact is assumed to occur on a single point on each surface, consistent to convex surfaces in contact. The pair of contact points is obtained by finding the minimum distance between these surfaces. For that, these points must satisfy the orthogonality conditions given by
$$\varvec{r}=\left[\begin{array}{c}{{\Gamma }}_{A {, }_{{\zeta }_{A}}}\bullet \varvec{g}\\ {{\Gamma }}_{A {, }_{{\theta }_{A}}}\bullet \varvec{g}\\ {{\Gamma }}_{A {, }_{{\zeta }_{B}}}\bullet \varvec{g}\\ {{\Gamma }}_{A {, }_{{\theta }_{B}}}\bullet \varvec{g}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right]$$
9
,
where comma “,” describes a partial derivative and \(\varvec{g}=\left({{\Gamma }}_{A}-{{\Gamma }}_{B}\right)\) is the gap function between the surfaces. It is worth noting that the gap function is dependent on the convective coordinates of each surface, but also on the surfaces’ vector position.
This set of nonlinear equations, the Local Contact Problem (LCP), is then solved by applying the Newton-Raphson Method or an optimization method. If penetration is verified, a pair of contact normal forces is imposed by a penalty method or a Lagrange multiplier. For more complex surfaces including concave regions, the best approach is to divide the surface into smaller convex surfaces as possible. Then the process consists of applying the same method on each surfaces pairs candidate to contact. A global search procedure can be used to first find which pairs are plausible to contact.
However, the model described in [8] takes advantage of a particular solution of this problem. Since the wheels are idealized as revolution surfaces and the rails as extrusion surfaces, one could restrict the analysis to the transversal plane of the track, and so to the contact between the 2D geometric profiles instead of the 3D complete surfaces, fixing one convective coordinate of each surface. Furthermore, if the geometric profiles are described by a set of circular arcs (Fig. 7), the LCP has an analytical solution, substantially reducing the problem’s complexity. This solution is based on the fact that the minimum distance between two circular arcs always occurs on the segment linking their centers. This distance can be easily calculated by the difference between the arc centers distance and the sum of arc radii.
Since this work analyzes the interaction of real wheelsets and tracks, the wheels and rails profiles are obtained by measurement instruments like MiniProf©, which provides a list of almost equidistant points describing the geometric profile. To transform these measurements into a set of tangent circular arcs (arc spline), Silva and Gay Neto [20] developed a fitting procedure that defines arcs based on a curvature function estimated from the points and a least-square optimization to fulfill the geometric and tolerance requirements. Thus, the two 2D geometric profiles on each wheel/rail pair are divided in multiple arcs candidates to contact. Then, multiple contact points can be obtained from the penetration among pairs of arcs, located on wheel and rail sides.
Despite this simplification, the model is still nonlinear due to the multiple contact points detection. The model using the dynamic relaxation method previously described is solved using the software GIRAFFE (Generic Interface Readily Accessible for Finite Elements) [21]. Then the results are post processed using the Hertz Contact Theory [22] to obtain the contact ellipse and pressure distribution, using as input the radii obtained by the arc spline fitting procedure.
Pummeling analysis
A pummelling analysis represents the accumulated effect of a set of consecutive contacts, such as profile wear or rolling contact fatigue. The pummelling analyses in this work are presented in terms of accumulated contact pressure. That is considered an indicator of the occurrence of rolling contact fatigue. Each pummelling analysis shows a graph of the summation of all individual contact pressure distributions over the rail profile. To plot this graph, each rail profile is linearized in a new coordinate S representing the length over the rail profile, as shown in Fig. 8. The coordinate S starts at zero at the profile field edge.
Four scenarios were simulated: tangent track with new rails, curved track with new rails, tangent track with worn rails, and curved track with worn rails. For each scenario, the contact of ten unique wheelsets were evaluated. Therefore, each pummelling analysis represents the summation of ten distinct contact pressure distributions. The wheels and rail profiles were obtained from the measurement of actual wheels and rails, the same used by Higa, Kina and Gay Neto [8] (Fig. 9 and Fig. 10). It includes wheel profiles with different characteristics, including thin flange, hollow, false flange, and no defects.
This study uses a section of a Brazilian heavy haul railway that runs 905 km from two state capitals in the country, Belo Horizonte and Vitória. The line has a metric gauge and mixed traffic includes both passenger and freight trains. Iron-ore is the main product transported.
For both SIMPACK® and VAMPIRE® simulations, tangent and curved track are evaluated in a single simulation. Thus, 20 simulations were performed in each software to evaluate all scenarios. The wagon travels at a constant speed of 65 km/h (18.05 m/s) in a track comprising of a tangent section with 150 m, a transition spiral for entry into the curve with 61.7 m, a curved section to the left with 195.6 m, a transition with 61.7 m and a final tangent section of 150 m; totaling 619 m traveled (Fig. 11). The track gauge is 1.0 m with rail cant of 1:40. At the curved section, the radius is 371.6 m and the superelevation is 57.2 mm. The total mass of the loaded wagon is 110,000 kg and the back-to-back distance between the wheels is 0.915 m.
For these simulations, the pummeling analyses were carried out with the contact data extracted from the front wheelset at the end of the first tangent section and the end of the curved section. Therefore, it represents the steady state of each section. Since both SIMPACK® and VAMPIRE® do not provide the contact pressure distribution as an output, such distributions were admitted as following the Hertz Contact Theory [22], i.e., a semi-elliptical distribution (Fig. 12). After the dynamic simulation, maximum contact pressure, the semi-axes of the ellipse and position of the contact on the rail were exported to MATLAB®, and the contact pressure distribution was calculated according to Eq. (10). In this equation, p0 is the maximum contact pressure and b represents the semi-axis of the contact ellipse in the transverse direction of the track. P is the contact load.
$$p\left(s\right)={p}_{0}\sqrt{1-{\left(\frac{s}{b}\right)}^{2}}$$
10
$${p}_{0}=\frac{3P}{2\pi ab}$$
11
For the quasi-static model, each scenario is simulated individually and 40 simulation were performed to evaluate all scenarios. The track and wagon parameters were the same as described for the multibody simulations.