Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 3, 39, 2001 CONTACT, MOTION AND WEAR IN RAILWAY MECHANICS Kurt Frischmuth FB Mathematik, Universität Rostock, Germany e-mail: kurt@sun2.math.uni-rostock.de One of the most intriguing problems in railway mechanics is the mo- deling of contact between a rail and wheel. On the one hand, in the contact zone the friction forces are determined – which are essential for the behavior of railway vehicles as multibody systems, see Kalker (1990), True (1993). On the other hand, the power extended by those frictional forces on the surfaces in contact leads to abrasion of particles andthus towearphenomena, seeBrommundt(1996),Frischmuth(1996), Langemann (1999). The aimof this paper is to present recent results on the couplingbetween motion of a railwayvehicle and evolutionof contact surfaces due towear. We concentrateonawheelwith initiallyperfect geometry,which changes due to frictional power dissipation during rolling on a perfect track. Key words: vehicle dynamics, friction, wear, railwaymechanics, contact problems, multibody systems 1. Motion Typically, vehicles are described as multibody systems – we include embs (ElasticMultiBodySystems) in this notion.Positions of allmasseswhich con- stitute the system under consideration are defined in terms of a finite number of coordinates y (or sometimes p = p(t)). Those coordinates may be (pre- ferably) independent, in that case their number – the dimension of y – is referred to as degrees of freedom. Otherwise constraints may be imposed, e.g. in the case of so-called kinematically closed loops. In the first case the temporal changes of y are governed by ordinary differential equations of second order m ( y(t) ) ÿ(t)= f ( t,y(t), ẏ(t),u ( t,y(t), ẏ(t) ) ) (1.1) 508 K.Frischmuth where the temporal derivatives are denoted by the dot. The mass matrix m represents the mass distribution within the described system, the function f comprises all (generalized) forces acting on the system. The function u is introduced to allow for a control motion, e.g. imposing a desired speed on the center of mass. A constrained system can be transformed, in principle, into such a setting using Lagrange formalism. However, from a numerical point of view this is rather not desirable. A direct approach to a differential-algebraic system m ( y(t) ) ÿ(t)= f ( t,y(t), ẏ(t),λ,u ( t,y(t), ẏ(t) ) ) +∇g⊤λ (1.2) 0= g ( t,y(t),u ( t,y(t), ẏ(t) ) ) leads to more efficient algorithms, see Gear et al. (1985), Hairer and Wanner (1991), Simeon et al. (1991), Arnold and Frischmuth (1998). For a rigid railway wheel, for instance, we have six degrees of freedom. However, assuming that the wheel is in continuous contact with a rigid rail, the vertical shift (elevation) of the mass center can be skipped. Similarly, for a rigid wheelset with bothwheels in contact with rails the roll angle φ can be evaluated for a given lateral shift y and yaw angle ψ. In both cases rather complicated geometrical calculations are required to obtain the redundant position variables from the degrees of freedom.To speed up the calculations, especially for simulations in real time, one prefers to do those calculations off-line during preprocessing, see Arnold and Frischmuth (1998). Of course, this works best if the bodies involved have symmetries and do not change with time. If the elasticity of the bodies in contact cannot be neglected then the constraints are not adequate, in which case we have equations of motion of the form (1.1). If we allow for temporary separation of the contact couple we may have to switch between different descriptions. The above amounts to the fact that the modeling of motion of a railway vehicle is quite a complicated task. There exists a number of successful simula- tion packages where the mentioned problems have been appropriately simpli- fied and solved. It is worthwhile tomention that depending on themotivation some rather simple models proved most interesting, see e.g. Kaas-Petersen (1986), True (1993). For the purpose of studyingwear phenomena it is obvious that we need to be able to cover relatively long time intervals. This excludes too complicated models. On the other hand, we need someminimum spatial resolution, thus a certain level of complication cannot be avoided. Contact, motion and wear in railway mechanics 509 2. Contact Most important for the mbs models introduced in the previous section is certainly the definition of the force function f, and in the case of constraints also of the function g. For ideal constraints, in the case of the formulation (1.2)1, we have constraint forces in the normal direction to the constraint manifold (of unknown intensity, determined byLagrange multipliers λ). In the presence of friction tangential forces depend on normal forces, i.e. on λ. Moreover, the forces depend on such details of the contact geometry as the curvatures of the surfaces in the contact zone. This dependence is complicated enough for reasonably smooth surfaceswhich have contact in exactly one spot. In railwaymechanics, unfortunately, wehave to facemultiple contact andhigh derivatives. Typically, thedetermination of contact forces is carried out in several steps. First, the geometrical contact point is determined. This is done in a purely geometric way, given the actual geometry of the contact partners and their positions in terms of mbs coordinates. For unworn profiles this dependence is well studied and available e.g. in the form of tables. Second, the so-called normal problem is solved. 3. Wear For this paper we assume that, given the multibody position vector y(t) and the generalized velocities ẏ(t), we are able to calculate the distributions of wear-relevant factors over the contact surfaces. The most frequently used wear model assumes the removed mass to be proportional to the dissipated energy. All we need in this case is local slip and tangential stresses in the contact region. Their scalar product, multiplied by the wear coefficient β, is the local speed of surface retreat. Direct numerical time integration of this quantity, for each point of the involved surfaces, in order to obtain the current positions of the worn sur- faces, fails for the obvious reason that integration time is too long and the integrand too oscillating. A frequently used technique is to amplify the real wear coefficient (by several orders), and shorten the integration time in the same way. This approach was criticized by Langemann (1999). Here we fol- low a different route.We postulate an evolution problem for the wear surface 510 K.Frischmuth under consideration of the general form ẋ(t,σ)=F(t,σ,x(t, ·))~n(t,σ) (3.1) where σ – surface parameter ~n – outer normal and −F is the wear speed. The latter may depend on time and position on the surface, but the most essential is the dependence on the actual surface as a whole. In the previous work we have made efforts to express this dependence, at least approximately, in an analytic way. Simulations with a nonlinear regu- larized inverse diffusion equation have led to instabilities of the wear process similar to the observed behavior. In this paper, however, we want to evaluate F by numerical simulation of amultibody system described by (1.1). Roughly speaking, the evaluation of F consists in calculating the integral mean values over time of the local wear for a mbs with frozen geometry. It turns out that – for the systems we studied – such mean values converged to unique limits, depending on the control motion (e.g. speed) and actual geometry x(t, ·), but independentof the initial values assumed for the degrees of freedom y in (1.1). This unique limit is assumed to be the postulatedwear speed −F ( t, ·,x(t, ·) ) . Efforts to give analytical representations of such amapping x(t, ·)→ ẋ n (t, ·) (3.2) acting essentially from a function space defined on the surface to itself, re- mained restricted to very simple setups, cf Brommundt (1996), Langemann (1999). Numerical evaluation of F depends ’only’ on our ability to describe the geometry and forces for worn surfaces represented by x(t, ·), to integrate equations of motion (1.1) and to integrate along the factors determining the wear. The remaining part is now the integration of evolution equation (3.1). Doing this by Euler’s explicit method is much like the naive algorithms based on amplification factors. This factor is the ratio between simulation time and step size of Euler’smethod.However, given formulation (3.1) we have now the choice of the whole wealth of known integration methods. In particular, we get a much better control of artificial oscillations. Contact, motion and wear in railway mechanics 511 4. Numerical study In the case of railway mechanics typically force and, if any, constraint functionsaregiven in the formofhundredsof lines of a source code.Thismakes analytical considerations extremely difficult, the more so as useful properties are not to be expected. In this paper our main interest is focused on the numerical integration of wear equation (3.1).Todo thiswe choose adynamicalmodel that is reasonably realistic, but on the other hand simple enough to allow quick and reliable integration of the auxiliary problems. We base this single wheelset model on the work by True (1993) but introduce some modifications in order to get feedback from the geometrical changes in the profiles. For the wear part, we choose the simplest possible law, assuming the loss of mass proportional to the dissipated frictional energy. We review now briefly the important ingredients of the dynamical model, stressing the modifications that come with time-dependent profiles. To start with, let us have a look at the standard rail and wheel profiles (UIC 60/S1002), cf Fig.1, as they are being met in Europe on new tracks, respectively vehicles. Each of the profiles is defined piecewise, second derivati- ves suffer jumps where the pieces meet. The rail profile is concave, the wheel profile has changing signs of curvature.Most important, the distance between the rail andwheel, which depends on their relative displacement, is essentially non-convex. In general, it has lots of local minima. Gravity causes, as long as it dominates the dynamical forces, theminimumof the distance function to be zero, the contact forces are applied at the spot where this minimum is assu- med. Hence, obviously, the distance function, and its minimizer are essential constituents of the problem formulation (1.1), respectively (1.2). Fig. 1. Rail and wheel profiles 512 K.Frischmuth By design, for all reasonable lateral and angular displacements, there sho- uld be a reasonably large region where the rail and wheel profiles come close to each other, cf Fig.1. Otherwise the stresses would exceed the yield stress. Zooming in shows, cf Fig.2, that for chosen position coordinates, there is a unique minimizer at −2cm from the origin in the lateral wheel surface di- rection. This point x c is called the point of geometrical contact. Under load, however, there is an elastic approach between the rail and wheel of the order 10−4m. The bottom of the graph of the actual distance function is flat for a patch of an extension of the order 10−2m. Obviously, the geometric contact point x c will not be exactly in the center of such a patch. Fig. 2. Distance function We have to stress that the point of contact and related quantities can be expressed in terms of the surface parameters σ of the wheel or railhead, or in coordinates of the 3Dphysical space. In the plotswe indicate by the subscripts r, w or ws whether we refer to the rail, wheel or wheelset. The meaning of x and y depends on the context, e.g. y2 is the generalized mbs coordinate of the rolling wheelset, y r is the lateral parameter of the railhead frame. The point of geometrical contact is determined by the generalized coor- dinates defining the position of the wheel in the system of the railhead. For symmetry reasons the lateral displacement, roll and yaw angles alone define Contact, motion and wear in railway mechanics 513 the contact point and the vertical displacement (compliance neglected), cf Ar- nold and Netter (1998). For frozen yaw angle the dependence is depicted in Fig.3. Note the step regions in the plot which correspond to the jumps of the contact point, e.g. from the tread to flange. This very unpleasant effect is absent in the case of conical wheels. While numerically attractive, especially for the dae-approach, this case is far too unrealistic for a serious study of we- ar. In fact, the lack of the flange leads to large lateral displacements so that the dissipated energy would spread over a region much wider than the actual wheel is. −0.02 −0.01 0 0.01 0.02 −3 −2 −1 0 1 2 3 −0.06 −0.04 −0.02 0 0.02 0.04 y w [m] φ [°] σ c [m ] Fig. 3. Point of geometrical contact Note that for each yaw angle ψ we get a different dependence.Moreover, like the location of the globalminimumthis characteristic depends in a sensiti- veway on the underlying profiles. Usually, these geometric data are calculated during preprocessing and interpolated at running time. For wear calculations, however, the procedure has to be repeated at each time step of the integration of (3.1). For a wheelset with a rigid axle the permanent contact to rails, see Fig.4, introduces a constraint, which eliminates the roll angle φ from the equations. In fact, shifting the wheelset to the left or right lifts the corresponding wheel 514 K.Frischmuth Fig. 4.Wheelset on rails −0.02 −0.01 0 0.01 0.02 0 1 2 3 −1 −0.5 0 0.5 1 y ws [m] ψ [°] φ [° ] Fig. 5. Roll angle of wheelset Contact, motion and wear in railway mechanics 515 −0.02 −0.01 0 0.01 0.02 0 0.5 1 1.5 2 2.5 3 −0.01 0 0.01 0.02 y ws [m] ψ [°] z w s [m ] Fig. 6. Vertical displacement up. Thuswe get a dependence of the roll angle (Fig.5) and the center of mass (Fig.6) on the lateral shift y and yaw angle ψ of the rigid wheelset. Multiplied by the axle load, this gives potential of the conservative forces acting on the wheelset. We have a stabilizing effect like for the motion in a half-pipe. However, the width of this half-pipe is finite, there is no extension of this potential beyond the interval [−0.07,0.06]. For proper exploitation the trajectories should mostly stay within the flat bottom of the slide, cf Meinke (2000). Themost delicate part of the dynamicalmodel, however, is the calculation of frictional forces. In terms ofmultibody system coordinates we have the load and rigid body slip as the input quantities. It turns out that locally, within the contact patch, the slip velocity is very inhomogeneous. The patch splits into a stick and slip region the size of which is essential for the calculation of the resultant tangential force, which is needed for time-integration of the mbs-coordinates. This interference of a finer scale of resolution is usually tre- ated with a level of accuracy dictated by the computational goals. There are reasonable algebraic approximations (Vermeulen Johnson), one step discreti- zedmethods (Fastsim) or full numerical solutions (Contact), cf Kalker (1990), Kik and piotrowski (1996). Our choice follows themethod byTrue (1993), but with the input variables taken from the actual worn profiles instead of the frozen values of the original model. 516 K.Frischmuth −0.02 −0.01 0 0.01 0.02 −0.02 −0.01 0 0.01 0.02 0 5 10 15 x 10 8 y r [m]xr [m] p [ N /m 2 ] Fig. 7. Normal stress in contact region −0.01 −0.005 0 0.005 0.01 −0.01 −0.005 0 0.005 0.01 y r [m] x r [ m ] Fig. 8. Tangential stress in contact region Contact, motion and wear in railway mechanics 517 In the first step, we calculate the distribution of the normal stress p (pres- sure), see Fig.7, followingHertz’s theory, see Hertz (1982). Thuswe disregard the tangential stresses and the known fact that the distance function is assu- med to be quadratic in Hertz’s theory. Obviously, a quadratic approximation is not valid in a domain big enough for the given order of normal loads and hence elastic approaches, cf Fig.2 again. However, better results come at a much higher numerical cost, so for our present purpose we accept the error in the normal component. In the consequent step, on the basis of the obtained normal stress and without any feedback, the local slip velocity and tangential stresses are appro- ximated. A sample result is presented in Fig.8. Only in a part of the patch we have slip, and thus, full contribution of µp to the tangential forces q. For a very high rigid body slip the whole patch becomes the slip region, then we have saturation and the friction force becomes µP with load L. For a moderate slip we get just a fraction of that limit value. For quick calculations a spline approximation to the full solution has to suffice. Now we have – at least in a very abbreviated way – described the major contributions to the wheelset model being sensitive to geometric changes of the wheel profiles. Let us assume given initial conditions for all free position variables and their time-derivatives. As the control motion we assume a pre- scribed longitudinal position of the center of mass (or of the end of a spring pulling the center of mass) of the wheelset. Then we have a well-posed initial value problem for the system (1.1). For a reasonable load, prescribed travel speed and initial conditions the solution exists for all the time (stays in the half-pipe). Stationary solutions are unstable above the so-called critical speed, solutions with flange contact are shown in Fig.9. andFig.10, cf Frischmuth et al. (1996). For certain speeds chaosmay be observed, cfKaas-Petersen (1986). In terms ofmultibodyvariables, given those solutions, we can calculate the dissipated power. In order to study wear, we need again the higher resolution of the local distribution within the contact patch. The scalar product of the local slip velocity and tangential stress gives the local surface intensity of the power dissipation. Obviously, this is a highly concentrated and quickly varying quantity. For- tunately, it turns out that taking an integral mean value over integration time tends very fast to a constant limit distribution, under given constant condi- tions. Within this paper, we assume that the wheel stays round. For out of round wheels we refer the Reader to works by Frischmuth (1996, 1997), Fri- schmuth and Langemann (1998). In this case all we need is the curve along the profile giving the average dissipated power per unit length and unit time, 518 K.Frischmuth Fig. 9. Hunting wheelset Fig. 10. Yaw angle versus lateral displacement Contact, motion and wear in railway mechanics 519 cf Chudzikiewicz (2000), Frischmuth (2000). Multiplied by the constant wear factor β this is equivalent to the prescribed normal speed −F at which the profile is abrased. −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0 0.2 0.4 0.6 0.8 1 0 0.005 0.01 0.015 0.02 t/t f [1] σ [m] −F Fig. 11. Power of dissipated energy Now, eventually, we have full description of the numerical evaluation of the wear speed −F, i.e. the operator on the right-hand side of evolution equation (3.1) for the wheel geometry. The integration time of (3.1) turns out to be sensitive, even in a 1Dcase, to the profile, having assumed symmetry of thewheel. Usually, explicit difference methods tend to oscillations which are easy to understand: abrasion of the material in the contact region leads to increasing of the distance function, hence the contact region is going tomove away. Too large time steps together with too fine resolution of the profile geometry leads to alternating phases of the abrasion in the neighboring regions. For stable and fast methods for solving (3.1) we refer to Sethian (1996). Fig.11 presents the power dissipation over the lifetime of the wheel. The change in the wheel profile is not really visible, it is still of the order of the elastic approach.Nonetheless, given the steep dependence of the contact point on theprofiledata, this leads to clear variation of the contact conditionsduring thewear process. In particular, due towear of initially exposed spots, thewear eventually reaches places that are not initially in contact whatever the mbs 520 K.Frischmuth positions. The actual wear can be approximated from the data presented in Fig.11by integrating along the timeaxis,multiplyingwith the constant β and modifying the profile in the normal direction to the reference configuration. 5. Conclusions Railway mechanics is a very wide field of research, the problem of rail- wheel contact is just one of many of its aspects, but certainly one of themost complicated. It cannot be uncoupled from other sub-problems, like models of the track, the subgrade, the whole vehicle. In this paper we assumed fargoing simplifications in order to concentrate on the numerical strategy for solving the surface evolution. Other important effects remain unaccounted for, e.g. plastic deformations, damage below the surface, vibrations, cf Bogacz and Dżula (1993). Ourmain point is that wear should bemodeled as an evolution problem, and appropriate numericalmethods should be used to calculate time dependence of the surfaces exposed to abrasive wear. References 1. ArnoldM., FrischmuthK., 1998, Solving Problems with Unilateral Con- straints byDAEMethods,Mathematics and Computers in Simulation, 47, 47- 67 2. Frischmuth K., Arnold M.. Hänler M., Netter H., 1996, Differential- gleichungen und singuläre Mannigfaltigkeiten in der dynamischen Simulation von Rad-Schiene-Systemen, in: Hoffmann, Jäger, Lohmann, Schunck (Edit.), Mathematik – Schlüsseltechnologie für die Zukunft, Proceedings des BMBF- Statusseminars München 1995, Springer-Verlag Berlin/Heidelberg 3. Arnold M., Netter H., 1998, Approximation of Contact Geometry in the Dynamical Simulation of Wheel-Rail Systems, Mathematical and Computer Modelling of Dynamical Systems, 4, 162-184 4. Bogacz R., Dżula S., 1993,Dynamics and Stability of aWheelset inRolling Contact Motion on Rails, in: Proc. of ITTG International Symposium on the Technological Innovation in Guided Transports, 871-883, Lille 5. Brommundt E., 1996, A SimpleMechanism for the Polygonalization of Rail- wayWheels byWear,Mechanics Research Communications Basic and Applied Contact, motion and wear in railway mechanics 521 6. Chudzikiewicz A., 2000, Zużycie kół kolejowych, in: Proceedings of the 7th Conference on Computer Simulation in Research andDevelopment, Kościelisko 7. Gear C.W., Gupta G.K., Leimkuhler B., 1985, Automatic Integration of Euler-Lagrange Equations with Constraints, J. Comp. Appl. Math., 12&13, 77-90 8. Hairer E., Wanner G., 1991, Solving Ordinary Differential Equations. II. Stiff and Differential-Algebraic Problems, Springer-Verlag, Berlin, Heidelberg, NewYork 9. Hertz H., 1882, Über die Berührung fester elastischer Körper, J. Reine und Angewandte Mathematik, 92 10. Kalker J.J., 1990, Three-Dimensional Elastic Bodies in Rolling Contact, Vol.2 of Solid Mechanics and its Applications, Kluwer Academic Publisher, Dordrecht 11. FrischmuthK., 1996,On theNumerical Solution ofRail-WheelContactPro- blems, J. Theoretical and Applied Mechanics, 34, 1, 7-15 12. FrischmuthK., 1997,Rail-WheelContact,Dynamical Simulation andDama- ge, In: W. Kosiński, R. de Boer, D. Gross (Edit.), Problems of Environmental and Damage Mechanics, IPPTPAN,Warszawa, 149-156 13. Frischmuth K., Langemann D., 1998, Numerical Analysis of Long-Term WearModels,Machine Dynamics Problems, 20, 113-122 14. Frischmuth K., 2000, On the Temporal Evolution of Wear Profiles, In: Z. Strzyżakowski andP.Lesiak (Edit.),Proceedings of the 4th TransCompCon- ference, Zakopane 15. Kaas-Petersen Ch., 1986, Chaos in a Railway Bogie, Acta Mechanica, 61, 89-107 16. Kik W., Piotrowski J., 1996, A Fast, Approximate Method to Calculate Normal Load at Contact Between Wheel and Rail and Creep Forces During Rolling, in: I. Zobory (Edit.)Proc. of the 2nd MiniConference on Contact Me- chanics and Wear of Rail/Wheel Systems, 52-61, TUBudapest 17. Langemann D., 1999,Numerische Analyse abrasiv verschleißender Mechani- scher Systeme, Fortschritt-Berichte 12.392, VDI-Verlag Düsseldorf 18. Meinke P., 2000, Über das Entgleisen von Radsätzen, In: K. Frischmuth (Edit.), Proceedings of the 7th Workshop on Dynamics of Rail/Wheel Dyna- mics, Rostock 19. SimeonB., FührerC., Rentrop P., 1991,Differential-Algebraic Equations in Vehicle SystemDynamics, Surveys on Mathematics for Industry, 1, 1-37 20. Sethian J.A., 1996, Level Set Methods, Cambridge University Press, Cam- bridge 522 K.Frischmuth 21. True H., 1993, Dynamics of a Rolling Wheelset, Applied Mechanics Reviews, 46, 7 Kontakt, ruch i zużycie w mechanice kolejowej Streszczenie Jednym z najbardziej ciekawych zagadnień mechaniki kolejowej jest modelowa- nie kontaktu między kołem a szyną. Z jednej strony w strefie kontaktu określone są siły tarcia – które są istotne dla zachowania pojazdów kolejowych jako układów wielu ciał, patrz Kalker (1990), True (1993). Z drugiej strony zaś moc dyssypowana przez siły tarcia na powierzchniach kontaktowych prowadzi do usuwania cząsteczek i tym samym do zjawisk zniszczenia, patrz Brommundt (1996), Frischmuth (1996), Langemann (1999). Celem niniejszej pracy jest prezentacja aktualnych wyników na temat sprzężenia między ruchempojazdukolejowegoa ewolucjąpowierzchni kontaktuwskutek zużycia. Dlaustalenia uwagi rozpatrujemykoło opoczątkowo idealnej geometrii, która zmienia się pod wpływem mocy sił tarcia dyssypowanej podczas toczenia po perfekcyjnym torze. Manuscript received February 21, 2001; accepted for print March 20, 2001