- Email: [email protected]
Finite element analysis of the frictional wheel-rail rolling contact using explicit and implicit methods
Author’s Accepted Manuscript Finite element analysis of the frictional wheel-rail rolling contact using explicit and implicit methods Moncef Toumi, Hugues Chollet, Honoré Yin
www.elsevier.com/locate/wear
PII: DOI: Reference:
S0043-1648(16)30122-3 http://dx.doi.org/10.1016/j.wear.2016.06.008 WEA101715
To appear in: Wear Received date: 6 October 2015 Revised date: 1 June 2016 Accepted date: 8 June 2016 Cite this article as: Moncef Toumi, Hugues Chollet and Honoré Yin, Finite element analysis of the frictional wheel-rail rolling contact using explicit and implicit methods, Wear, http://dx.doi.org/10.1016/j.wear.2016.06.008 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Finite element analysis of the frictional wheel-rail rolling contact using explicit and implicit methods Moncef Toumia,b,c,*, Hugues Cholleta, Honoré Yinb a b
Department of COSYS/GRETTIA, IFSTTAR, 14-20 Boulevard Newton, Champs sur Marne, France Laboratoire Navier, UMR8205, École des Ponts, IFSTTAR, CNRS, UPE, Champs-sur-Marne, France c IRT Railenium, F-59300, Famars, France *[email protected]
ABSTRACT A three-dimensional finite element rolling contact model between wheel and rail with increased spin effect is developed in this paper to study the normal and the tangential contact problems, as well as the creep force characteristics in elasticity and in elasto-plasticity. Three finite element analyses using the explicit and the implicit integration scheme are employed and the finite element solution is compared to the solution of the CONTACT software in elasticity. The elastic solution with the implicit methods is sensitive to the energy dissipation related to the time increment size. The optimum time increment size is found as a compromise between computational efficiency and calculation accuracy. The numerical results show a good agreement between the finite element methods and CONTACT for the normal and the tangential contact solution. The initial slopes of the creep force characteristics in the longitudinal and the lateral direction from the implicit solution are slightly lower than the results from CONTACT. When the elasto-plastic material behavior is applied, the maximum contact pressure and the von Mises stresses are reduced for both explicit and implicit methods and the plastic deformation occurs in all the contact area. The contact patch becomes egg shaped instead of the elliptic and the area of stick zone area increases. In addition, the initial slopes of the creepage-creep force curves are reduced by 15% in the longitudinal direction and by 13% in the lateral direction. Keywords Rolling contact, rail-wheel tribology, friction, finite element modelling, creep-force curve.
1. Introduction The contact area between the wheel and the rail is subjected to high pressure and shear stresses that could directly lead to their damage. The calculation of the contact forces acting over this contact area is also an important step in the determination of the dynamic performance of vehicle-track system. Several authors such as Carter [1], Johnson [2,3,4] and Kalker [5,6,7] studied these contact forces for a quasi-static movement of the axle and following their investigations, many rolling contact theories and software tools are now available to calculate the tangential forces in the contact area, expressing them as a function of the creepage which is the relative velocity between the wheel and the rail and the famous coefficients of stiffness provided theoretically by Kalker’s table [5]. These coefficients were measured on test-benches several times between the 20’s and 80’s for a longitudinal, lateral and spin creepage [8-11]. In 1967, Hobbs [9] made a survey of most of these experimental works showing that some of those measurements present a decrease in the creepage-creep force curves that can reach 50% compared to Kalker’s theoretical prediction. This dispersion is confirmed few years later by Nayak et al. [10] in 1970 who studied the influence of some factors on friction and creep in rolling contact such as rolling velocity, surface contamination, surface roughness and dynamic load due to irregular track. Surface contamination is mentioned to be the most important factor when it comes to the large difference between experimental measurements and Kalker theoretical curves. Then, this dispersion is also confirmed by Pearce and Rose [11] in 1982 who presented measurements of the creepage-creep force relationship made on service tracks. The effect of elasto-plasticity is generally disregarded by different authors. Several authors developed models that are able to describe the reduced creep coefficients and the falling friction for higher creepage values. Vollebregt [12] has modeled the reduction in the initial slope by adding an interfacial layer between the wheel and the rail. He implemented his work in the software CONTACT based on Kalker’s complete theory [7] using the boundary element method (BEM). Meierhofer et al. [13] developed a model based on the half-space assumption to describe the influence of a Third Body Layer on the traction characteristic in the longitudinal direction. Many authors used the finite element (FE) method to model the wheel-rail rolling contact [14-16]. Most of the works
conducted so far aim essentially to perform an elasto-plastic stress analysis. However, FE studies that analyze the creep force characteristics and the parameters that could influence its initial slopes and the Kalker’s creep coefficients seem to lack in the literature. In this paper, the FE method is used to model the frictional rolling contact between the wheel and the rail under elastic and elasto-plastic conditions. The study focuses on the solution concerning the contact stresses and the creep force characteristics in the longitudinal and the lateral direction. The result in elasticity is compared with the Hertz theory and CONTACT [17]. The effect of elasto-plasticity on the contact shear stresses and on the initial slope of the creep force curves is investigated. Different FE analyses (explicit dynamic analysis, implicit moderate energy dissipation analysis and implicit quasi-static analysis) are conducted and compared to each other in this work.
2. Numerical modeling The commercial software ABAQUS [18] is used in this work to model the rolling contact between the wheel and the rail in elasticity then in elasto-plasticity. The Newmark- method is used in this code to solve the general differential equations and to find the dynamic response. Depending on and some other parameters values, the time integration scheme could be chosen. The explicit and the implicit integration schemes were used in this work. The explicit solver is very effective in the case of high nonlinearities but the solution stability depends on the time step. In fact, the time step required to guarantee the stability with explicit time integration must satisfy the Courant stability condition [19]. This condition requires that the time step is small enough to ensure that a sound wave may not cross the smallest element during one time step. This time step is proportional to the smallest mesh size in the model. Thus, it requires a very small time step to reach the stability since the minimum mesh size in our model is sometimes under 0.5 mm which leads to a high calculation time. The implicit integration scheme has the advantage of being unconditionally stable for linear systems. However, the solution stability for nonlinear problem as the rolling wheel-rail contact problem might be sensitive to the numerical energy dissipation related to the size of the time increment. As long as the size of the time increment is small enough, the energy dissipation should be sufficient to avoid some fluctuations that could affect the normal and the tangential contact solution. For a nonlinear problem, the computational cost for explicit integration scheme rises linearly with the reduction of the time increment size while it rises more rapidly than linearly when the implicit integration scheme is used. Therefore, to reduce the computational cost and meet good accuracy with respect to creep force characteristics, one should choose the appropriate FE approach between implicit and explicit schemes. The implicit dynamic analysis is performed for two applications: the «Moderate Dissipation (MD)» and the «Quasi-Static (QS)» application. These two procedures possess some form of numerical dissipation that allows damping of higher vibrations modes without affecting too strongly the lower modes. The MD procedure uses an extension of the Newmarkmethod, the Hilber-Hughes-Taylor time integration scheme. This algorithm allows a better control of the numerical dissipation by parameters other than the time step [20]. The value of these parameters set in ABAQUS for this application [18] allows a moderate amount of numerical dissipation. The QS procedure uses the backward Euler time integration scheme [18]. This scheme tends to be more dissipative than the Hilber-Hughes-Taylor scheme. A considerable numerical dissipation may be used to obtain convergence while the computational cost is minimized by taking large time increments when possible. The numerical dissipation allowed by these two implicit procedures could be combined by others forms of energy dissipation like plasticity or viscous damping.
Fig. 1. FE mesh of the wheel-rail (S1002-UIC60) rolling contact model. In this work, the phenomena occurring close to the contact area are treated. The upper part of the wheel and the bottom part of the rail has been removed from the 3D model in order to reduce the total number of elements and thus reduce the computational cost. The wheel and the rail have been meshed with 8-node hexahedron elements. The mesh is refined throughout the tread (Fig. 1) defined by the two planes around the contact point which is found by an analytic procedure available on the wheel-rail calculation platform Expro [21]. Zhao et al. [15] showed that a 1.3×1.3 mm mesh size can provide a satisfactory accuracy for most engineering applications. Thus, for the sake of accuracy, the smallest mesh size for which detailed contact results were studied is 0.5 mm with a total element number of 380000 while the smallest one for which the creep force curves were drawn is 1 mm with a total element number of 91000. The wheel profile in Fig. 1 is a conical S1002 wheel type (the tread was replaced by a cone) with a taper of 1/10 in order to increase the spin effect, and the rail is an UIC60 type tilted by 1/10. The length of the rail is 250 mm and the length of the wheel is 260 mm. The two surfaces in contact are considered as smooth. The relative sliding between the wheel and the rail is defined using the «finite sliding» formulation. The «surface-to-surface» algorithm is used to calculate the interaction between the wheel and the rail in the presence of a coefficient of friction =0.2. The penalty friction formulation is employed in this paper1. The material of both wheel and rail is assumed to have material properties according to the R900A steel commonly used for the manufacture of railway rails with a Young's modulus = 210 GPa and a Poisson's ratio ν= 0.28, the density of steel is = 7800 kg/m3. All parameters of the model are shown in Table 1. Table 1 : Parameters of the model Parameters Values 7.5 kN Load per wheel 50 km/h Forward speed 210 GPa Young’s modulus 0.28 Poisson ratio 7800 kg/m3 Density 489 MPa Yield stress 0.2 Coefficient of friction 920 mm Rolling wheel diameter Rail head radius 300 mm Wheel conicity 1/10 Tilt of the rail 1/10 Rail length 250 mm The value of some parameters in this work like lateral position of the wheel with respect to the rail , tilt of the rail and wheel conicity might not exist in a revenue traffic wheel-rail contact situation. They have been selected only to ensure a centered elliptic contact area (according to Hertzian conditions) so that the results could easily be compared with other methods. A Cartesian coordinate system ( ) is defined as shown in Fig. 1, the center O is connected to the axle center, the wheel is tied to the rigid axle (not represented in Fig. 1) where the normal load is applied at its center as indicated. The general analysis was made in two steps: during step 1 the normal contact was established and the normal load was applied. The rail rests on a rigid foundation and the wheel, in adhesion with the axle, moves in the direction where the load is applied as reported in Fig. 1. The step 2 is the rolling step: A forward speed was applied to the rail in the opposite sense of rolling. The rolling speed was applied around the axis of rotation of the axle in the sense indicated in Fig. 1. This axis coincides with the axis for a yaw angle equal to zero. The speeds and were applied gradually using the « Smooth step » procedure [18] and reach their maximum values after the first 13.5 % of the traveled distance d3= 195 mm, see Fig. 7. Both values of and were set in order to control the value of the longitudinal creepage given for rigid bodies by: (1) where is the rolling radius that depends on the lateral position . During this same step, the rotation of the wheel around is released and the advance of the rail along axis begins from a distance d1=0 mm. Thus, the wheel and the rail can respectively roll and move forward. The normal load is maintained 1
The Lagrange multiplier method has been used in some cases with the same results but with slightly higher computational cost
until the end of the simulation. The lateral creepage is defined by imposing a yaw angle α around the axis i.e. rotation angle between the wheel and the rail, see Fig. 1. Due to this angle, the rolling speed of the wheel will have two components and respectively along the and the axis, which are given by: (2) (3) The spin creepage is handled automatically by the calculation procedure due to the geometry of bodies in contact i.e. the conicity of the wheel and the curvature of the rail and due to the rail inclination. Unless otherwise specified, the detailed contact solution is investigated at the distance d=115 mm for = -0.7‰ and for =0 mrad. The FE model is further developed in the present work by considering an elasto-plastic behavior of the R900A material.
3. Instabilities associated with the FE implicit method As stated in the previous section, the implicit integration schemes are unconditionally stable for linear systems. However, the study of the wheel-rail rolling contact includes different sources of nonlinearities which are mainly due to frictional contact. In this case, some irregularities or instabilities related to the high frequency vibrations may arise and could affect the solution accuracy. Automatic time incrementation is used by default in the software ABAQUS for nonlinear implicit procedures. The adjustment of the time increment size is controlled by two factors: the convergence behavior of the Newton-Raphson iterations and the accuracy of the time integration. Using the implicit MD and QS procedures to study the wheel-rail rolling contact in elasticity leads to constant time increment =3.7e-04 s.
Fig. 2. FE solutions of the von Mises stresses with the implicit MD procedure during rolling; on top: at contact rail surface, at the bottom: in depth (longitudinal section view of the rail, deformation scale amplified by factor 100); (a) time increment =3.7e-04 s, (b) =3.7e-04 s and Rayleigh damping, (c) =1e-04 s. Fig. 2(a) shows the von Mises stresses with during the rolling step when the implicit DM procedure is used. The top view shows an asymmetry of the von Mises stresses distribution relative to the plane (yoz) with a stress concentration at the rear of the contact surface. This asymmetry can also be observed in the longitudinal section view of the contact pressure in Fig. 3 that also becomes more elongated in the longitudinal direction.
Fig. 3. Longitudinal section view of the contact pressure from the FE implicit MD solution. In order to better understand the origin of this asymmetry, the longitudinal section view of the rail is presented in Fig. 2(a) (cutting plane A-A) and the deformation was amplified by a factor 100. The result shows waves on the rear of the contact area with an approximately wavelength of 10 mm. The most likely origin of these waves is a dynamic response in the structure that arises under the effect of tangential traction when the dynamic equations are solved using the implicit MD procedure. The moderate energy dissipation associated with the Hilber-Hughes-Taylor scheme seems not enough to damp the oscillation at the rear of the contact surface that leads to asymmetry in the pressure and von Mises stresses distribution (Fig. 2(a)). To improve the damping behavior with the implicit MD method, two techniques are used. First, in addition to the artificial damping, the Rayleigh damping [22] with coefficients and is introduced in the model. The Rayleigh constants are used to set the amount of damping. Their values =0.01 and =1e-06 are evaluated from an empirical procedure after testing different values. Fig. 2(b) presents the von Mises stresses when the Rayleigh coefficients are applied. The longitudinal section view (cutting plane B-B) shows that the oscillation is completely damped and the waves at the rear of the contact area disappear which give a symmetrical von Mises stresses and pressure distribution (Fig. 3) as expected by theory [4]. The second technique consists on reducing the time increment size, the amount of energy dissipation related to the time step should then increase. Fig. 2(c) shows the von Mises stresses with a time increment about 4 times smaller than . As with Rayleigh damping, the waves at the rear of the contact area (cutting plane C-C) disappear which also gives a symmetrical distribution for the von Mises stresses. Fig. 3 shows that the contact pressure with is in good agreement with the Rayleigh damping solution. It is also shown from Fig. 2 that the additional damping with these two techniques reduces the maximum von Mises stresses that decreases from =1002 MPa with the time increment to =777 MPa when the Rayleigh damping is applied and =801 MPa when the time increment is used. It should be mentioned that these oscillations don’t exist with the implicit QS procedure. This may be explained by the amount of the numerical dissipation related to the backward Euler scheme that is naturally higher than that with the implicit MD procedure. A second form of irregularities related to a different mode of vibration is observed in the contact shear stresses. The 3D view in Fig. 4(d) shows that these contact shear stresses oscillate inside the contact surface when the default time increment is used. Neither the Rayleigh damping nor the implicit QS procedure had succeeded in damping these oscillations. In contrast, the 3D view of the contact shear stress on the right of the Fig. 4 shows that by reducing the time increment size, the amplitude and the frequency of these oscillations decrease until completely disappear with =2e-05 s. This can also clearly be seen on the left of the Fig. 4 for the normalized contact shear stress with respect to , with is the normal nodal pressure.
Fig. 4. Left: normalized contact shear stress with respect to (b),(e) , (c),(f) . = -0.7‰, =0 mrad. =0.2. Table 2 : Relative difference with respect to Maximum pressure Semi-axis a Semi-axis b Maximum von Mises stresses Maximum contact shear stress Longitudinal resultant force
-2.5% 10.3% 0% 31.5% 8% 27.2%
( /
); Right: 3D view of contact shear stress; (a),(d)
,
for different quantities 0.3% 0% 0% 3.7% 5% 7%
Reducing the time increment size leads to a higher calculation time which is about 29.5 h with , 10.5 h with and 7.5 h with for eight 2.8 GHz processors. Table 2 shows the difference between results with , and . The difference is significant between results with and while the difference is small between results with and . So, will be considered as the optimal time increment to meet both a good accuracy and reasonable calculation time. It will then be used in this paper to study the creep curves characteristics. However, to ensure a better comparison with the reference solution given by CONTACT for purpose of validation, the oscillation in Fig. 4 should be minimized. Thus, the detailed normal and tangential contact results for the FE implicit methods are presented in the following section with the time increment .
4. Results and discussion A 3D-FE model of frictional rolling contact between a conical S1002 wheel type and an UIC60 rail type has been developed using the procedure described above. In the following, the results related to the normal and tangential problem will be presented first under elastic conditions then under elasto-plastic conditions. 4.1. Elastic solution 4.1.1. Normal contact solution Fig. 5 shows the distribution of contact pressure during the rolling step with different FE solutions and with CONTACT. The contact patch from FE solutions is elliptic as we expected. The maximum pressure and the semi-axis a and b are in good agreement with Hertz theory and with CONTACT as shown in Table 3.
Fig. 5. Distribution of normal contact pressure given by different numerical methods; (a) CONTACT, (b) Explicit dynamic, (c) Implicit QS, (d) Implicit MD. Table 3 : Normal contact solution with different methods
CONTACT Hertz Explicit Dyn Implicit QS Implicit MD
Semi-axis a (mm)
Semi-axis b (mm)
6.75 6.66 7.25 7.00 7.25
5.25 5.05 5.50 5.50 5.50
Maximum pressure (MPa) 1132 1160 1126 1158 1158
The previous comparison shows that all FE methods give a very satisfactory accuracy with respect to CONTACT for the normal solution. However, the computational cost is clearly higher when the implicit QS and MD methods are used with the time increment which last approximately 27.5 h compared to 9.5 h with the explicit dynamic method. Indeed, unlike the explicit method, the implicit time integration schemes need to inverse the operator matrix at each time increment which is computationally expensive. Thus, when only the normal solution is of interest, it is suggested to use the explicit method or the implicit method with a larger time increment size such as in order to save unnecessary computational cost. 4.1.2. Tangential contact solution Fig. 6 shows a comparison between the FE solutions from different methods and CONTACT for the normalized contact shear stress with respect to (left) and for the contact shear stress contour and vectors (right). The perfect agreement between the implicit QS and MD solutions in Fig. 6 is due to the same time increment used with these two methods, which applies a same amount of artificial damping. The vectors of contact shear stress in Fig. 6(right) shows the trend of the torque due to spin to rotate the wheel around the rail relative to the normal axis . The contact shear stress dissymmetry according to the longitudinal axis (y=0) is due to this spin effect. Both explicit and implicit methods are in agreement with CONTACT. The vectors are oriented in the clockwise sense. The contour of contact shear stress with the FE implicit methods shows a distribution close to that given by CONTACT with a similar concentration at the rear of the contact area. However, the FE explicit method shows a larger zone where the highest contact shear stress contours are concentrated as shown in Fig. 6(f). This difference has no effect on the maximum contact shear stress that is approximately equal to =210 MPa with all numerical methods. This result is important for wheel-rail contact modelling since the development of wear and the crack initiation path are determined by the magnitude of the contact shear stress. The stick and slip regions are deduced from the normalized contact shear stress with respect to according to Coulomb law (the bodies stick when < and slip when = ). The longitudinal creepage combined with the spin creepage leads to the final stick/slip area presented in Fig. 6(left). Again, the FE implicit solutions are in good agreement with CONTACT. However, the contact shear stress with the explicit dynamic solution saturates in larger zone than CONTACT. The slip region is at the leading edge as shown in Fig. 6(b). The normal and tangential contact solution showed that the FE implicit QS and MD methods give practically the same results since the time increment is small enough to apply the same numerical dissipation. Thus, in the following sections, the implicit method will be represented only by the implicit MD procedure in order to avoid further unnecessary computational cost. The study of the creep force characteristics with the FE method is computationally expensive. In choosing the appropriate FE approach, the best compromise between the computational cost and the solution accuracy has to be considered. Based on
the discussions above and the time increment sensitivity conducted in section 3, it can be concluded that the FE implicit method with time increment =1e-04 s gives the best choice. It will then be used in the following section to study the creep force characteristics in elasticity and in elasto-plasticity. The FE explicit method will be represented only in the longitudinal direction in order to check the consistency between different numerical solutions.
Fig. 6. Left: normalized contact shear stress with respect to ( / ); Right: contour and vectors of contact shear stress; (a),(e) CONTACT, (b),(f) Explicit dynamic, (c),(g) Implicit QS, (d),(h) Implicit MD. =-0.7‰, =0 mrad. =0.2. 4.1.3. Creep force characteristics Accurate estimation of the creep force characteristics is of extreme importance for a vehicle dynamic simulation. The FE model presented in section 2 is used to study these creep force characteristics in elasticity and in elasto-plasticity. The creep force and are calculated by summing the nodal contact shear force respectively in the longitudinal and the lateral direction. Fig. 7 represents these longitudinal and lateral creep forces with respect to time from the FE explicit and implicit MD solutions during the rolling step. The longitudinal creepage value =-0.7‰ leads to the negative longitudinal force which corresponds to a braking operational condition. The steady state is rapidly reached with the two methods when the rolling velocity reaches her maximum value at t2=0.003 s after rolling over the distance d2=41 mm. It is then maintained until the end of the simulation at the distance d3=195 mm.
Fig. 7. Longitudinal and lateral creep force with respect to time with the explicit and the implicit MD methods. =0 mrad. =0.2. Lower figure: position along the rail.
= -0.7‰,
The expression of the longitudinal creepage reported in eq.(1) is based on the assumption of rigidity of the two bodies in contact. The local deformations of the wheel and the rail bodies due to their elasticity are not taken into account in this expression. When the contact between the wheel and the rail take place, the wheel undergoes an elastic deformation and the rolling radius , initially calculated with the rigid-body assumption and used in CONTACT, decreases slightly from to . Therefore, the FE creepage value is different from the one applied in CONTACT and an offset on the creepage value will be observed on the creep force characteristics between the FE solutions and CONTACT. In order to facilitate the comparison between the different FE solutions and CONTACT, the origins of different longitudinal creep force curves in Fig. 8(a) and Fig. 12(a) are set at the same value (0,0). Thus, curves from the FE solutions are shifted by their offset values. The creep force characteristics in the longitudinal and the lateral direction from the FE methods are represented in Fig. 8 and compared with CONTACT. A slight difference is observed between the implicit solution and CONTACT in the longitudinal direction (Fig. 8(a)). The saturation lower limit is in agreement with =-3‰ and the saturation upper limit with the implicit method is slightly higher with =4‰ compared to =3.5‰ with CONTACT. The initial slope of this curve with the FE implicit solution is slightly lower than CONTACT. The highest initial slope in the longitudinal direction is given by the explicit FE solution where the lower and upper saturation levels are reached for respectively =-2‰ and =2.5‰.
Fig. 8. The creep force characteristics in elasticity from the FE methods and comparison with CONTACT; (a) longitudinal creep force curve (b) lateral creep force curve. Fig. 8(b) shows that the initial slope of the creep force curve in the lateral direction is again slightly lower with the FE implicit method than with CONTACT. The saturation lower limit with the implicit method and CONTACT is respectively =-6 mrad and =-5 mrad and the saturation upper limit is respectively =3.5 mrad and =2 mrad. According to Kalker [5], the lateral force is not only due to the lateral creepage but also to the spin as illustrated in the following equation: (
)
(4)
where G is the shear modulus and are the Kalker’s coefficients. The spin can therefore lead to a lateral force even if the yaw angle is zero. This fact is verified in the FE solutions. The presence of a relatively high lateral force =5053 N in Fig. 8(b) for =0 mrad is explained by this spin effect. Multibody dynamic codes use generally the Kalker’s coefficients Cij to calculate directly the contact tangential forces using Kalker’s linear theory. From this study, more accurate coefficients might be estimated from the initial slope of the creep force curves given by the FE solution in Fig. 8. In the following section, the effect of the elasto-plastic material behavior on the normal and tangent contact solution and on the creep force characteristics is studied. 4.2. Elasto-plastic solution One of the advantages of the FE method is its capability to handle nonlinear material behavior like elasto-plastic material. The normal wheel-rail contact problem has been treated in recent years in plasticity. For example, Wiest et al. [23] used the FE method to calculate the contact pressure, the contact patch size and the penetration depth under elasto-plastic conditions. Zhao et al. [14] solved the normal and the tangential problems using a FE explicit transient analysis. Besides Zhao’s solution and Wen et al. work [16], FE studies for frictional rolling contact between wheel and rail in elasto-plasticity seems still to lack in literature. In this section, the wheel-rail contact during one rolling cycle is studied in elasto-plasticity. The normal and the tangential contact solution, as well as the creep force characteristics are studied. A nonlinear material model with isotropic hardening is used in this study. The material data from the R900A steel are available from a tensile test. The yield stress is =489 MPa. The program CONTACT could not be used in this study since it is limited to contact between elastic bodies. The comparison will then be made only between the different FE solutions. 4.2.1. von Mises stresses at the surface Fig. 9 shows the FE solutions for the von Mises stresses at the rail contact surface during the rolling step and under elastoplastic conditions. The result shows a good agreement between the explicit and the implicit FE solutions. According to this criterion and considering the yield stress of the R900A steel, the plastic flow occurs on nearly all the contact area. The maximum von Mises stresses is reduced to 491 MPa with both FE methods. Moreover, residual stress is observed in Fig. 9
at the rail surface where plastic deformation occurs after the wheel passage. It should be mentioned that additional energy dissipation due to plastic deformation is applied under elasto-plastic conditions. This dissipation contributes to damp the von Mises stresses and the contact shear stress fluctuations observed respectively in Fig. 2(a) and Fig. 4(a).
Fig. 9. Von Mises stresses with elasto-plastic material behavior from: (a) Explicit dynamic solution, (c) Implicit MD solution. 4.2.2. Normal contact solution Under the effect of plastic deformation at the contact surface, the deformed wheel-rail profiles are slightly modified at the contact level compared to the initial undeformed profiles. Thus, the contact pressure distribution from different FE solutions is no longer symmetric relative to (yoz) plane, see Fig. 10. The contact patch changes from elliptic shape under the elastic conditions (Fig. 5) to egg shape. The maximum pressure is reduced to 987 MPa and 994 MPa respectively from explicit dynamic solution and the implicit MD solution. The normal solutions from the FE explicit and implicit methods under elasto-plastic conditions are in good agreement.
Fig. 10. Distribution of normal contact pressure in elasto-plasticity: (a) Explicit dynamic solution, (b) Implicit MD solution. 4.2.3. Tangential contact solution Fig. 11 shows the the normalized contact shear stress with respect to ( / ) (left) and the contact shear stress contour and vectors (right) obtained from different FE solutions in elasto-plasticity. The solutions obtained by the implicit and the explicit dynamic methods (Fig. 11(a),(b)) show a larger stick area compared to the elastic solution (Fig. 6(b),(d)). This can be explained by the wheel plastic deformation. When the plastic deformation at the wheel contact surface occurs, the rolling radius of the wheel is reduced slightly by . Here, is due to the elastic deformation and is due to the plastic deformation at the wheel contact surface. This additional reduction leads to a rolling radius that is smaller than the rolling radius presented in the elastic solution (see section 4.1.3). Actually, the longitudinal creepage calculated from eq.(1) considers the bodies in contact as rigid. However, when the FE method is used, the wheel and the rail undergo elasto-plastic deformation and eq.(1) becomes inaccurate. This contributes in the difference observed first between the FE elastic solutions and the FE elasto-plastic solutions in Fig. 11 and second between the FE elastic solutions and CONTACT solution in Fig. 6 since the latter uses the same equation but don’t take into account the change in rolling radius under the elastic deformation. Comparing the elasto-plastic solution in Fig. 11(c),(d) to the elastic one in Fig. 6(f),(h), it is seen that elasto-plasticity has a negligible effect on the vectors of contact shear stress. However, the maximum contact shear stress is slightly reduced from =210 MPa to =191 MPa.
Fig. 11. Left: normalized contact shear stress with respect to ( / ); Right: contour and vectors of contact shear stress; (a),(c) Explicit dynamic solution, (b),(d) Implicit MD solution. =-0.7‰, =0 mrad. =0.2. 4.2.4. Creep force characteristics The effect of elasto-plasticity on the creep force characteristics in the longitudinal and the lateral direction is studied using the FE implicit MD analysis. Fig. 12(a) represents the longitudinal force with respect to the creepage . The curve shows that the elasto-plastic behavior of the steel reduces the initial slope of the creep force curve by 15 % and consequently reduces the Kalker’s coefficient. However, this same behavior has no effect on the saturation region of the curve since the saturation zone depends mainly on the coefficient of friction and not on the material behavior. Similarly, Fig. 12(b) shows a reduction by 13 % in the initial slope of the curve representing the lateral force with respect to the creepage under the elasto-plastic behavior which reduces the Kalker’s coefficient.
Fig. 12. The creep force characteristics in elasto-plasticity from the implicit method and comparison with the elastic solution; (a) longitudinal creep force curve, (b) lateral creep force curve.
5. Conclusions A three dimensional FE model has been developed to study the frictional wheel-rail rolling contact with increased spin effect in elasticity and elasto-plasticity with three different FE analyses: the explicit dynamic analysis, the implicit quasistatic analysis and the implicit moderate energy dissipation analysis. The solution related to the normal and the tangential contact problems, as well as to the creep force characteristics is studied. The validation in elasticity is made essentially by comparing different FE solutions and CONTACT. The conclusions from this study are as follows: (1) The normal contact solution in elasticity from different FE methods is in good agreement with Hertz theory and CONTACT for the maximum contact pressure, the contact patch size and the contact shape. The tangential contact solution in elasticity from the FE implicit methods is in good agreement with CONTACT for the contour and vectors of the contact shear stress as well as for the stick/slip area. The FE explicit dynamic solution shows a slip area in the leading edge of the contact surface. However, the value of the maximum contact shear stress is not affected by the FE method used and agrees well with CONTACT. (2) Under the elasto-plasticity effect, the maximum von Mises stresses and the maximum pressure are reduced and the shape of the contact patch changes from elliptic to egg shaped. Furthermore, the stick area becomes larger with both explicit and implicit methods and the maximum contact shear stress is slightly reduced. (3) The initial slope of the creep force characteristics obtained from the FE implicit solution under elastic conditions is slightly lower than CONTACT solution in the longitudinal and the lateral direction. Under elasto-plastic conditions, it is shown for a realistic example that the initial slopes of these curves are reduced by 15% in the longitudinal direction and by 13% in the lateral direction compared to the FE elastic solution. However, the relevance of these values should
be further investigated in more situations such as different wheel-rail couples in order to decide which reduction coefficients to use in vehicle dynamics software. (4) The FE implicit solution is sensitive to the energy dissipation related to the size of the time increment. Indeed, if the time increment is not small enough, oscillations at the rear of the contact area as well as inside the contact shear stress are observed which leads to stick/slip distribution widely different compared to CONTACT. Thus, the time increment must be carefully chosen by analyzing time increment sensitivity at the beginning of the study. Acknowledgments
The authors thank the Technological Research Institute RAILENIUM for their funding through the project CERVIFER and TATA STEEL for providing the material data References [1] F.W. Carter, On the action of a locomotive driving wheel, Proc. R. Soc. Lond. 112 (1926) 151-157. [2] K.L. Johnson, The effect of a tangential contact force upon the rolling motion of an elastic sphere on a plane, J. Appl. Mech., Trans ASME, 25 (1958) 338-340. [3] P.J. Vermeulen and K. L. Johnson, Contact of non-spherical bodies transmitting tangential forces, J. Appl. Mech., 31 (1964) 338-340. [4] K.L. Johnson, Contact Mechanics, Cambridge University Press, Cambridge, 1987. [5] J. J. Kalker, On the rolling contact of two elastic bodies in presence of dry friction, (Ph.D. thesis), Delft University of Technology, 1967. [6] J.J. Kalker, A fast algorithm for the simplified theory of rolling contact, Vehicle Syst. Dyn., 11 (1982) 1-13. [7] J.J. Kalker, Three-Dimensional Elastic Bodies in Rolling Contact, Kluwer Academic Publishers, Dordrecht, 1990. [8] C. Pritchard, Traction between rolling steel Surfaces - A Survey of Railway and Laboratory Experiments. Proceedings of the 7th Leeds/Lyon Symposium on Tribology, 1980. [9] A. Hobbs, A survey of creep, British Railways Research Dpt Report note DYN 52, 1967. [10] P.R. Nayak, S. Hariharan, R. Stern, R. Abilock and P.A. March, Friction and creep in rolling contact, Federal Railroad Administration, Dpt of transportation, 1970. [11] T.G. Pearce and K.A. Rose, Tangential force-creepage relationships in theory and practice, Contact mechanics and wear of rail/wheel systems, University of Waterloo Press, 1982. [12] E.A.H. Vollebregt, Numerical modeling of measured railway creep versus creep-force curves with CONTACT, Wear 314 (2014) 87-95. [13] A. Meierhofer, C. Hardwick, R. Lewis, K. Six and P. Dietmaier, Third body layer-experimental results and a model describing its influence on the traction coefficient, Wear 314 (2014) 148-154. [14] X. Zhao, Z. Li and R. Dollevoet, Solution of the wheel rail rolling contact in elasticity and elasto-plasticity using a transient finite element, International Conference on Contact Mechanics and Wear of Rail/Wheel Systems (CM2009), Firenze, ltaly, 2009. [15] X. Zhao and Z. Li, The solution of frictional wheel-rail rolling contact with a 3D transient finite element model: Validation and error analysis, Wear, 271 (2011) 444-452. [16] Z. Wen, L. Wu, W. Li, X. Jin and M. Zhu, Three-dimensional elastic–plastic stress analysis of wheel–rail rolling contact, Wear, 271 (2011) 426-436. [17] E.A.H. Vollebregt, User Guide for CONTACT, Vollebregt & Kalker's Rolling and Sliding Contact Model, Version13.1, VORtech, 2013. [18] Abaqus User Guide, Version 6.12, Dassault Systèmes Simulia Corp., 2012. [19] R. Courant, K. Friedrichs and H. Lewy, On the partial differential equations of mathematical physics, Math. Ann. 100 (1928) 32-74. [20] H.M. Hilber, T.J.R. Hughes and R.L. Taylor, Improved numerical dissipation for time integration algorithms in structural dynamics, Earthquake engineering and structural dynamics, 5 (1977) 283-292. [21] P. Aknin, H. Chollet, M. Sebes, Expro, Exploration de profils roue-rail, Version 4.0, INRETS, 2005. [22] P. Jehel, Modélisation numérique des phénomènes d'amortissement par dissipation d'énergie matérielle dans les structures de type portique en béton armé sous séisme (Ph.D. thesis), Ecole Normale Supérieure de Cachan. 2009. [23] M. Wiest, E. Kassac, W. Davesa, J.C.O. Nielsen and H. Ossberger, Assessment of methods for calculating contact pressure in wheel-rail/switch contact, Wear 265 (2008) 1439-1445.
Highlights
A 3D FE model is developed to study the frictional rolling wheel-rail contact
Explicit and implicit integration schemes are used and the solutions are compared
Implicit solution is sensitive to energy dissipation related to time increment size
The elasto-plasticity reduces the initial slopes of the creep force characteristics
www.elsevier.com/locate/wear
PII: DOI: Reference:
S0043-1648(16)30122-3 http://dx.doi.org/10.1016/j.wear.2016.06.008 WEA101715
To appear in: Wear Received date: 6 October 2015 Revised date: 1 June 2016 Accepted date: 8 June 2016 Cite this article as: Moncef Toumi, Hugues Chollet and Honoré Yin, Finite element analysis of the frictional wheel-rail rolling contact using explicit and implicit methods, Wear, http://dx.doi.org/10.1016/j.wear.2016.06.008 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Finite element analysis of the frictional wheel-rail rolling contact using explicit and implicit methods Moncef Toumia,b,c,*, Hugues Cholleta, Honoré Yinb a b
Department of COSYS/GRETTIA, IFSTTAR, 14-20 Boulevard Newton, Champs sur Marne, France Laboratoire Navier, UMR8205, École des Ponts, IFSTTAR, CNRS, UPE, Champs-sur-Marne, France c IRT Railenium, F-59300, Famars, France *[email protected]
ABSTRACT A three-dimensional finite element rolling contact model between wheel and rail with increased spin effect is developed in this paper to study the normal and the tangential contact problems, as well as the creep force characteristics in elasticity and in elasto-plasticity. Three finite element analyses using the explicit and the implicit integration scheme are employed and the finite element solution is compared to the solution of the CONTACT software in elasticity. The elastic solution with the implicit methods is sensitive to the energy dissipation related to the time increment size. The optimum time increment size is found as a compromise between computational efficiency and calculation accuracy. The numerical results show a good agreement between the finite element methods and CONTACT for the normal and the tangential contact solution. The initial slopes of the creep force characteristics in the longitudinal and the lateral direction from the implicit solution are slightly lower than the results from CONTACT. When the elasto-plastic material behavior is applied, the maximum contact pressure and the von Mises stresses are reduced for both explicit and implicit methods and the plastic deformation occurs in all the contact area. The contact patch becomes egg shaped instead of the elliptic and the area of stick zone area increases. In addition, the initial slopes of the creepage-creep force curves are reduced by 15% in the longitudinal direction and by 13% in the lateral direction. Keywords Rolling contact, rail-wheel tribology, friction, finite element modelling, creep-force curve.
1. Introduction The contact area between the wheel and the rail is subjected to high pressure and shear stresses that could directly lead to their damage. The calculation of the contact forces acting over this contact area is also an important step in the determination of the dynamic performance of vehicle-track system. Several authors such as Carter [1], Johnson [2,3,4] and Kalker [5,6,7] studied these contact forces for a quasi-static movement of the axle and following their investigations, many rolling contact theories and software tools are now available to calculate the tangential forces in the contact area, expressing them as a function of the creepage which is the relative velocity between the wheel and the rail and the famous coefficients of stiffness provided theoretically by Kalker’s table [5]. These coefficients were measured on test-benches several times between the 20’s and 80’s for a longitudinal, lateral and spin creepage [8-11]. In 1967, Hobbs [9] made a survey of most of these experimental works showing that some of those measurements present a decrease in the creepage-creep force curves that can reach 50% compared to Kalker’s theoretical prediction. This dispersion is confirmed few years later by Nayak et al. [10] in 1970 who studied the influence of some factors on friction and creep in rolling contact such as rolling velocity, surface contamination, surface roughness and dynamic load due to irregular track. Surface contamination is mentioned to be the most important factor when it comes to the large difference between experimental measurements and Kalker theoretical curves. Then, this dispersion is also confirmed by Pearce and Rose [11] in 1982 who presented measurements of the creepage-creep force relationship made on service tracks. The effect of elasto-plasticity is generally disregarded by different authors. Several authors developed models that are able to describe the reduced creep coefficients and the falling friction for higher creepage values. Vollebregt [12] has modeled the reduction in the initial slope by adding an interfacial layer between the wheel and the rail. He implemented his work in the software CONTACT based on Kalker’s complete theory [7] using the boundary element method (BEM). Meierhofer et al. [13] developed a model based on the half-space assumption to describe the influence of a Third Body Layer on the traction characteristic in the longitudinal direction. Many authors used the finite element (FE) method to model the wheel-rail rolling contact [14-16]. Most of the works
conducted so far aim essentially to perform an elasto-plastic stress analysis. However, FE studies that analyze the creep force characteristics and the parameters that could influence its initial slopes and the Kalker’s creep coefficients seem to lack in the literature. In this paper, the FE method is used to model the frictional rolling contact between the wheel and the rail under elastic and elasto-plastic conditions. The study focuses on the solution concerning the contact stresses and the creep force characteristics in the longitudinal and the lateral direction. The result in elasticity is compared with the Hertz theory and CONTACT [17]. The effect of elasto-plasticity on the contact shear stresses and on the initial slope of the creep force curves is investigated. Different FE analyses (explicit dynamic analysis, implicit moderate energy dissipation analysis and implicit quasi-static analysis) are conducted and compared to each other in this work.
2. Numerical modeling The commercial software ABAQUS [18] is used in this work to model the rolling contact between the wheel and the rail in elasticity then in elasto-plasticity. The Newmark- method is used in this code to solve the general differential equations and to find the dynamic response. Depending on and some other parameters values, the time integration scheme could be chosen. The explicit and the implicit integration schemes were used in this work. The explicit solver is very effective in the case of high nonlinearities but the solution stability depends on the time step. In fact, the time step required to guarantee the stability with explicit time integration must satisfy the Courant stability condition [19]. This condition requires that the time step is small enough to ensure that a sound wave may not cross the smallest element during one time step. This time step is proportional to the smallest mesh size in the model. Thus, it requires a very small time step to reach the stability since the minimum mesh size in our model is sometimes under 0.5 mm which leads to a high calculation time. The implicit integration scheme has the advantage of being unconditionally stable for linear systems. However, the solution stability for nonlinear problem as the rolling wheel-rail contact problem might be sensitive to the numerical energy dissipation related to the size of the time increment. As long as the size of the time increment is small enough, the energy dissipation should be sufficient to avoid some fluctuations that could affect the normal and the tangential contact solution. For a nonlinear problem, the computational cost for explicit integration scheme rises linearly with the reduction of the time increment size while it rises more rapidly than linearly when the implicit integration scheme is used. Therefore, to reduce the computational cost and meet good accuracy with respect to creep force characteristics, one should choose the appropriate FE approach between implicit and explicit schemes. The implicit dynamic analysis is performed for two applications: the «Moderate Dissipation (MD)» and the «Quasi-Static (QS)» application. These two procedures possess some form of numerical dissipation that allows damping of higher vibrations modes without affecting too strongly the lower modes. The MD procedure uses an extension of the Newmarkmethod, the Hilber-Hughes-Taylor time integration scheme. This algorithm allows a better control of the numerical dissipation by parameters other than the time step [20]. The value of these parameters set in ABAQUS for this application [18] allows a moderate amount of numerical dissipation. The QS procedure uses the backward Euler time integration scheme [18]. This scheme tends to be more dissipative than the Hilber-Hughes-Taylor scheme. A considerable numerical dissipation may be used to obtain convergence while the computational cost is minimized by taking large time increments when possible. The numerical dissipation allowed by these two implicit procedures could be combined by others forms of energy dissipation like plasticity or viscous damping.
Fig. 1. FE mesh of the wheel-rail (S1002-UIC60) rolling contact model. In this work, the phenomena occurring close to the contact area are treated. The upper part of the wheel and the bottom part of the rail has been removed from the 3D model in order to reduce the total number of elements and thus reduce the computational cost. The wheel and the rail have been meshed with 8-node hexahedron elements. The mesh is refined throughout the tread (Fig. 1) defined by the two planes around the contact point which is found by an analytic procedure available on the wheel-rail calculation platform Expro [21]. Zhao et al. [15] showed that a 1.3×1.3 mm mesh size can provide a satisfactory accuracy for most engineering applications. Thus, for the sake of accuracy, the smallest mesh size for which detailed contact results were studied is 0.5 mm with a total element number of 380000 while the smallest one for which the creep force curves were drawn is 1 mm with a total element number of 91000. The wheel profile in Fig. 1 is a conical S1002 wheel type (the tread was replaced by a cone) with a taper of 1/10 in order to increase the spin effect, and the rail is an UIC60 type tilted by 1/10. The length of the rail is 250 mm and the length of the wheel is 260 mm. The two surfaces in contact are considered as smooth. The relative sliding between the wheel and the rail is defined using the «finite sliding» formulation. The «surface-to-surface» algorithm is used to calculate the interaction between the wheel and the rail in the presence of a coefficient of friction =0.2. The penalty friction formulation is employed in this paper1. The material of both wheel and rail is assumed to have material properties according to the R900A steel commonly used for the manufacture of railway rails with a Young's modulus = 210 GPa and a Poisson's ratio ν= 0.28, the density of steel is = 7800 kg/m3. All parameters of the model are shown in Table 1. Table 1 : Parameters of the model Parameters Values 7.5 kN Load per wheel 50 km/h Forward speed 210 GPa Young’s modulus 0.28 Poisson ratio 7800 kg/m3 Density 489 MPa Yield stress 0.2 Coefficient of friction 920 mm Rolling wheel diameter Rail head radius 300 mm Wheel conicity 1/10 Tilt of the rail 1/10 Rail length 250 mm The value of some parameters in this work like lateral position of the wheel with respect to the rail , tilt of the rail and wheel conicity might not exist in a revenue traffic wheel-rail contact situation. They have been selected only to ensure a centered elliptic contact area (according to Hertzian conditions) so that the results could easily be compared with other methods. A Cartesian coordinate system ( ) is defined as shown in Fig. 1, the center O is connected to the axle center, the wheel is tied to the rigid axle (not represented in Fig. 1) where the normal load is applied at its center as indicated. The general analysis was made in two steps: during step 1 the normal contact was established and the normal load was applied. The rail rests on a rigid foundation and the wheel, in adhesion with the axle, moves in the direction where the load is applied as reported in Fig. 1. The step 2 is the rolling step: A forward speed was applied to the rail in the opposite sense of rolling. The rolling speed was applied around the axis of rotation of the axle in the sense indicated in Fig. 1. This axis coincides with the axis for a yaw angle equal to zero. The speeds and were applied gradually using the « Smooth step » procedure [18] and reach their maximum values after the first 13.5 % of the traveled distance d3= 195 mm, see Fig. 7. Both values of and were set in order to control the value of the longitudinal creepage given for rigid bodies by: (1) where is the rolling radius that depends on the lateral position . During this same step, the rotation of the wheel around is released and the advance of the rail along axis begins from a distance d1=0 mm. Thus, the wheel and the rail can respectively roll and move forward. The normal load is maintained 1
The Lagrange multiplier method has been used in some cases with the same results but with slightly higher computational cost
until the end of the simulation. The lateral creepage is defined by imposing a yaw angle α around the axis i.e. rotation angle between the wheel and the rail, see Fig. 1. Due to this angle, the rolling speed of the wheel will have two components and respectively along the and the axis, which are given by: (2) (3) The spin creepage is handled automatically by the calculation procedure due to the geometry of bodies in contact i.e. the conicity of the wheel and the curvature of the rail and due to the rail inclination. Unless otherwise specified, the detailed contact solution is investigated at the distance d=115 mm for = -0.7‰ and for =0 mrad. The FE model is further developed in the present work by considering an elasto-plastic behavior of the R900A material.
3. Instabilities associated with the FE implicit method As stated in the previous section, the implicit integration schemes are unconditionally stable for linear systems. However, the study of the wheel-rail rolling contact includes different sources of nonlinearities which are mainly due to frictional contact. In this case, some irregularities or instabilities related to the high frequency vibrations may arise and could affect the solution accuracy. Automatic time incrementation is used by default in the software ABAQUS for nonlinear implicit procedures. The adjustment of the time increment size is controlled by two factors: the convergence behavior of the Newton-Raphson iterations and the accuracy of the time integration. Using the implicit MD and QS procedures to study the wheel-rail rolling contact in elasticity leads to constant time increment =3.7e-04 s.
Fig. 2. FE solutions of the von Mises stresses with the implicit MD procedure during rolling; on top: at contact rail surface, at the bottom: in depth (longitudinal section view of the rail, deformation scale amplified by factor 100); (a) time increment =3.7e-04 s, (b) =3.7e-04 s and Rayleigh damping, (c) =1e-04 s. Fig. 2(a) shows the von Mises stresses with during the rolling step when the implicit DM procedure is used. The top view shows an asymmetry of the von Mises stresses distribution relative to the plane (yoz) with a stress concentration at the rear of the contact surface. This asymmetry can also be observed in the longitudinal section view of the contact pressure in Fig. 3 that also becomes more elongated in the longitudinal direction.
Fig. 3. Longitudinal section view of the contact pressure from the FE implicit MD solution. In order to better understand the origin of this asymmetry, the longitudinal section view of the rail is presented in Fig. 2(a) (cutting plane A-A) and the deformation was amplified by a factor 100. The result shows waves on the rear of the contact area with an approximately wavelength of 10 mm. The most likely origin of these waves is a dynamic response in the structure that arises under the effect of tangential traction when the dynamic equations are solved using the implicit MD procedure. The moderate energy dissipation associated with the Hilber-Hughes-Taylor scheme seems not enough to damp the oscillation at the rear of the contact surface that leads to asymmetry in the pressure and von Mises stresses distribution (Fig. 2(a)). To improve the damping behavior with the implicit MD method, two techniques are used. First, in addition to the artificial damping, the Rayleigh damping [22] with coefficients and is introduced in the model. The Rayleigh constants are used to set the amount of damping. Their values =0.01 and =1e-06 are evaluated from an empirical procedure after testing different values. Fig. 2(b) presents the von Mises stresses when the Rayleigh coefficients are applied. The longitudinal section view (cutting plane B-B) shows that the oscillation is completely damped and the waves at the rear of the contact area disappear which give a symmetrical von Mises stresses and pressure distribution (Fig. 3) as expected by theory [4]. The second technique consists on reducing the time increment size, the amount of energy dissipation related to the time step should then increase. Fig. 2(c) shows the von Mises stresses with a time increment about 4 times smaller than . As with Rayleigh damping, the waves at the rear of the contact area (cutting plane C-C) disappear which also gives a symmetrical distribution for the von Mises stresses. Fig. 3 shows that the contact pressure with is in good agreement with the Rayleigh damping solution. It is also shown from Fig. 2 that the additional damping with these two techniques reduces the maximum von Mises stresses that decreases from =1002 MPa with the time increment to =777 MPa when the Rayleigh damping is applied and =801 MPa when the time increment is used. It should be mentioned that these oscillations don’t exist with the implicit QS procedure. This may be explained by the amount of the numerical dissipation related to the backward Euler scheme that is naturally higher than that with the implicit MD procedure. A second form of irregularities related to a different mode of vibration is observed in the contact shear stresses. The 3D view in Fig. 4(d) shows that these contact shear stresses oscillate inside the contact surface when the default time increment is used. Neither the Rayleigh damping nor the implicit QS procedure had succeeded in damping these oscillations. In contrast, the 3D view of the contact shear stress on the right of the Fig. 4 shows that by reducing the time increment size, the amplitude and the frequency of these oscillations decrease until completely disappear with =2e-05 s. This can also clearly be seen on the left of the Fig. 4 for the normalized contact shear stress with respect to , with is the normal nodal pressure.
Fig. 4. Left: normalized contact shear stress with respect to (b),(e) , (c),(f) . = -0.7‰, =0 mrad. =0.2. Table 2 : Relative difference with respect to Maximum pressure Semi-axis a Semi-axis b Maximum von Mises stresses Maximum contact shear stress Longitudinal resultant force
-2.5% 10.3% 0% 31.5% 8% 27.2%
( /
); Right: 3D view of contact shear stress; (a),(d)
,
for different quantities 0.3% 0% 0% 3.7% 5% 7%
Reducing the time increment size leads to a higher calculation time which is about 29.5 h with , 10.5 h with and 7.5 h with for eight 2.8 GHz processors. Table 2 shows the difference between results with , and . The difference is significant between results with and while the difference is small between results with and . So, will be considered as the optimal time increment to meet both a good accuracy and reasonable calculation time. It will then be used in this paper to study the creep curves characteristics. However, to ensure a better comparison with the reference solution given by CONTACT for purpose of validation, the oscillation in Fig. 4 should be minimized. Thus, the detailed normal and tangential contact results for the FE implicit methods are presented in the following section with the time increment .
4. Results and discussion A 3D-FE model of frictional rolling contact between a conical S1002 wheel type and an UIC60 rail type has been developed using the procedure described above. In the following, the results related to the normal and tangential problem will be presented first under elastic conditions then under elasto-plastic conditions. 4.1. Elastic solution 4.1.1. Normal contact solution Fig. 5 shows the distribution of contact pressure during the rolling step with different FE solutions and with CONTACT. The contact patch from FE solutions is elliptic as we expected. The maximum pressure and the semi-axis a and b are in good agreement with Hertz theory and with CONTACT as shown in Table 3.
Fig. 5. Distribution of normal contact pressure given by different numerical methods; (a) CONTACT, (b) Explicit dynamic, (c) Implicit QS, (d) Implicit MD. Table 3 : Normal contact solution with different methods
CONTACT Hertz Explicit Dyn Implicit QS Implicit MD
Semi-axis a (mm)
Semi-axis b (mm)
6.75 6.66 7.25 7.00 7.25
5.25 5.05 5.50 5.50 5.50
Maximum pressure (MPa) 1132 1160 1126 1158 1158
The previous comparison shows that all FE methods give a very satisfactory accuracy with respect to CONTACT for the normal solution. However, the computational cost is clearly higher when the implicit QS and MD methods are used with the time increment which last approximately 27.5 h compared to 9.5 h with the explicit dynamic method. Indeed, unlike the explicit method, the implicit time integration schemes need to inverse the operator matrix at each time increment which is computationally expensive. Thus, when only the normal solution is of interest, it is suggested to use the explicit method or the implicit method with a larger time increment size such as in order to save unnecessary computational cost. 4.1.2. Tangential contact solution Fig. 6 shows a comparison between the FE solutions from different methods and CONTACT for the normalized contact shear stress with respect to (left) and for the contact shear stress contour and vectors (right). The perfect agreement between the implicit QS and MD solutions in Fig. 6 is due to the same time increment used with these two methods, which applies a same amount of artificial damping. The vectors of contact shear stress in Fig. 6(right) shows the trend of the torque due to spin to rotate the wheel around the rail relative to the normal axis . The contact shear stress dissymmetry according to the longitudinal axis (y=0) is due to this spin effect. Both explicit and implicit methods are in agreement with CONTACT. The vectors are oriented in the clockwise sense. The contour of contact shear stress with the FE implicit methods shows a distribution close to that given by CONTACT with a similar concentration at the rear of the contact area. However, the FE explicit method shows a larger zone where the highest contact shear stress contours are concentrated as shown in Fig. 6(f). This difference has no effect on the maximum contact shear stress that is approximately equal to =210 MPa with all numerical methods. This result is important for wheel-rail contact modelling since the development of wear and the crack initiation path are determined by the magnitude of the contact shear stress. The stick and slip regions are deduced from the normalized contact shear stress with respect to according to Coulomb law (the bodies stick when < and slip when = ). The longitudinal creepage combined with the spin creepage leads to the final stick/slip area presented in Fig. 6(left). Again, the FE implicit solutions are in good agreement with CONTACT. However, the contact shear stress with the explicit dynamic solution saturates in larger zone than CONTACT. The slip region is at the leading edge as shown in Fig. 6(b). The normal and tangential contact solution showed that the FE implicit QS and MD methods give practically the same results since the time increment is small enough to apply the same numerical dissipation. Thus, in the following sections, the implicit method will be represented only by the implicit MD procedure in order to avoid further unnecessary computational cost. The study of the creep force characteristics with the FE method is computationally expensive. In choosing the appropriate FE approach, the best compromise between the computational cost and the solution accuracy has to be considered. Based on
the discussions above and the time increment sensitivity conducted in section 3, it can be concluded that the FE implicit method with time increment =1e-04 s gives the best choice. It will then be used in the following section to study the creep force characteristics in elasticity and in elasto-plasticity. The FE explicit method will be represented only in the longitudinal direction in order to check the consistency between different numerical solutions.
Fig. 6. Left: normalized contact shear stress with respect to ( / ); Right: contour and vectors of contact shear stress; (a),(e) CONTACT, (b),(f) Explicit dynamic, (c),(g) Implicit QS, (d),(h) Implicit MD. =-0.7‰, =0 mrad. =0.2. 4.1.3. Creep force characteristics Accurate estimation of the creep force characteristics is of extreme importance for a vehicle dynamic simulation. The FE model presented in section 2 is used to study these creep force characteristics in elasticity and in elasto-plasticity. The creep force and are calculated by summing the nodal contact shear force respectively in the longitudinal and the lateral direction. Fig. 7 represents these longitudinal and lateral creep forces with respect to time from the FE explicit and implicit MD solutions during the rolling step. The longitudinal creepage value =-0.7‰ leads to the negative longitudinal force which corresponds to a braking operational condition. The steady state is rapidly reached with the two methods when the rolling velocity reaches her maximum value at t2=0.003 s after rolling over the distance d2=41 mm. It is then maintained until the end of the simulation at the distance d3=195 mm.
Fig. 7. Longitudinal and lateral creep force with respect to time with the explicit and the implicit MD methods. =0 mrad. =0.2. Lower figure: position along the rail.
= -0.7‰,
The expression of the longitudinal creepage reported in eq.(1) is based on the assumption of rigidity of the two bodies in contact. The local deformations of the wheel and the rail bodies due to their elasticity are not taken into account in this expression. When the contact between the wheel and the rail take place, the wheel undergoes an elastic deformation and the rolling radius , initially calculated with the rigid-body assumption and used in CONTACT, decreases slightly from to . Therefore, the FE creepage value is different from the one applied in CONTACT and an offset on the creepage value will be observed on the creep force characteristics between the FE solutions and CONTACT. In order to facilitate the comparison between the different FE solutions and CONTACT, the origins of different longitudinal creep force curves in Fig. 8(a) and Fig. 12(a) are set at the same value (0,0). Thus, curves from the FE solutions are shifted by their offset values. The creep force characteristics in the longitudinal and the lateral direction from the FE methods are represented in Fig. 8 and compared with CONTACT. A slight difference is observed between the implicit solution and CONTACT in the longitudinal direction (Fig. 8(a)). The saturation lower limit is in agreement with =-3‰ and the saturation upper limit with the implicit method is slightly higher with =4‰ compared to =3.5‰ with CONTACT. The initial slope of this curve with the FE implicit solution is slightly lower than CONTACT. The highest initial slope in the longitudinal direction is given by the explicit FE solution where the lower and upper saturation levels are reached for respectively =-2‰ and =2.5‰.
Fig. 8. The creep force characteristics in elasticity from the FE methods and comparison with CONTACT; (a) longitudinal creep force curve (b) lateral creep force curve. Fig. 8(b) shows that the initial slope of the creep force curve in the lateral direction is again slightly lower with the FE implicit method than with CONTACT. The saturation lower limit with the implicit method and CONTACT is respectively =-6 mrad and =-5 mrad and the saturation upper limit is respectively =3.5 mrad and =2 mrad. According to Kalker [5], the lateral force is not only due to the lateral creepage but also to the spin as illustrated in the following equation: (
)
(4)
where G is the shear modulus and are the Kalker’s coefficients. The spin can therefore lead to a lateral force even if the yaw angle is zero. This fact is verified in the FE solutions. The presence of a relatively high lateral force =5053 N in Fig. 8(b) for =0 mrad is explained by this spin effect. Multibody dynamic codes use generally the Kalker’s coefficients Cij to calculate directly the contact tangential forces using Kalker’s linear theory. From this study, more accurate coefficients might be estimated from the initial slope of the creep force curves given by the FE solution in Fig. 8. In the following section, the effect of the elasto-plastic material behavior on the normal and tangent contact solution and on the creep force characteristics is studied. 4.2. Elasto-plastic solution One of the advantages of the FE method is its capability to handle nonlinear material behavior like elasto-plastic material. The normal wheel-rail contact problem has been treated in recent years in plasticity. For example, Wiest et al. [23] used the FE method to calculate the contact pressure, the contact patch size and the penetration depth under elasto-plastic conditions. Zhao et al. [14] solved the normal and the tangential problems using a FE explicit transient analysis. Besides Zhao’s solution and Wen et al. work [16], FE studies for frictional rolling contact between wheel and rail in elasto-plasticity seems still to lack in literature. In this section, the wheel-rail contact during one rolling cycle is studied in elasto-plasticity. The normal and the tangential contact solution, as well as the creep force characteristics are studied. A nonlinear material model with isotropic hardening is used in this study. The material data from the R900A steel are available from a tensile test. The yield stress is =489 MPa. The program CONTACT could not be used in this study since it is limited to contact between elastic bodies. The comparison will then be made only between the different FE solutions. 4.2.1. von Mises stresses at the surface Fig. 9 shows the FE solutions for the von Mises stresses at the rail contact surface during the rolling step and under elastoplastic conditions. The result shows a good agreement between the explicit and the implicit FE solutions. According to this criterion and considering the yield stress of the R900A steel, the plastic flow occurs on nearly all the contact area. The maximum von Mises stresses is reduced to 491 MPa with both FE methods. Moreover, residual stress is observed in Fig. 9
at the rail surface where plastic deformation occurs after the wheel passage. It should be mentioned that additional energy dissipation due to plastic deformation is applied under elasto-plastic conditions. This dissipation contributes to damp the von Mises stresses and the contact shear stress fluctuations observed respectively in Fig. 2(a) and Fig. 4(a).
Fig. 9. Von Mises stresses with elasto-plastic material behavior from: (a) Explicit dynamic solution, (c) Implicit MD solution. 4.2.2. Normal contact solution Under the effect of plastic deformation at the contact surface, the deformed wheel-rail profiles are slightly modified at the contact level compared to the initial undeformed profiles. Thus, the contact pressure distribution from different FE solutions is no longer symmetric relative to (yoz) plane, see Fig. 10. The contact patch changes from elliptic shape under the elastic conditions (Fig. 5) to egg shape. The maximum pressure is reduced to 987 MPa and 994 MPa respectively from explicit dynamic solution and the implicit MD solution. The normal solutions from the FE explicit and implicit methods under elasto-plastic conditions are in good agreement.
Fig. 10. Distribution of normal contact pressure in elasto-plasticity: (a) Explicit dynamic solution, (b) Implicit MD solution. 4.2.3. Tangential contact solution Fig. 11 shows the the normalized contact shear stress with respect to ( / ) (left) and the contact shear stress contour and vectors (right) obtained from different FE solutions in elasto-plasticity. The solutions obtained by the implicit and the explicit dynamic methods (Fig. 11(a),(b)) show a larger stick area compared to the elastic solution (Fig. 6(b),(d)). This can be explained by the wheel plastic deformation. When the plastic deformation at the wheel contact surface occurs, the rolling radius of the wheel is reduced slightly by . Here, is due to the elastic deformation and is due to the plastic deformation at the wheel contact surface. This additional reduction leads to a rolling radius that is smaller than the rolling radius presented in the elastic solution (see section 4.1.3). Actually, the longitudinal creepage calculated from eq.(1) considers the bodies in contact as rigid. However, when the FE method is used, the wheel and the rail undergo elasto-plastic deformation and eq.(1) becomes inaccurate. This contributes in the difference observed first between the FE elastic solutions and the FE elasto-plastic solutions in Fig. 11 and second between the FE elastic solutions and CONTACT solution in Fig. 6 since the latter uses the same equation but don’t take into account the change in rolling radius under the elastic deformation. Comparing the elasto-plastic solution in Fig. 11(c),(d) to the elastic one in Fig. 6(f),(h), it is seen that elasto-plasticity has a negligible effect on the vectors of contact shear stress. However, the maximum contact shear stress is slightly reduced from =210 MPa to =191 MPa.
Fig. 11. Left: normalized contact shear stress with respect to ( / ); Right: contour and vectors of contact shear stress; (a),(c) Explicit dynamic solution, (b),(d) Implicit MD solution. =-0.7‰, =0 mrad. =0.2. 4.2.4. Creep force characteristics The effect of elasto-plasticity on the creep force characteristics in the longitudinal and the lateral direction is studied using the FE implicit MD analysis. Fig. 12(a) represents the longitudinal force with respect to the creepage . The curve shows that the elasto-plastic behavior of the steel reduces the initial slope of the creep force curve by 15 % and consequently reduces the Kalker’s coefficient. However, this same behavior has no effect on the saturation region of the curve since the saturation zone depends mainly on the coefficient of friction and not on the material behavior. Similarly, Fig. 12(b) shows a reduction by 13 % in the initial slope of the curve representing the lateral force with respect to the creepage under the elasto-plastic behavior which reduces the Kalker’s coefficient.
Fig. 12. The creep force characteristics in elasto-plasticity from the implicit method and comparison with the elastic solution; (a) longitudinal creep force curve, (b) lateral creep force curve.
5. Conclusions A three dimensional FE model has been developed to study the frictional wheel-rail rolling contact with increased spin effect in elasticity and elasto-plasticity with three different FE analyses: the explicit dynamic analysis, the implicit quasistatic analysis and the implicit moderate energy dissipation analysis. The solution related to the normal and the tangential contact problems, as well as to the creep force characteristics is studied. The validation in elasticity is made essentially by comparing different FE solutions and CONTACT. The conclusions from this study are as follows: (1) The normal contact solution in elasticity from different FE methods is in good agreement with Hertz theory and CONTACT for the maximum contact pressure, the contact patch size and the contact shape. The tangential contact solution in elasticity from the FE implicit methods is in good agreement with CONTACT for the contour and vectors of the contact shear stress as well as for the stick/slip area. The FE explicit dynamic solution shows a slip area in the leading edge of the contact surface. However, the value of the maximum contact shear stress is not affected by the FE method used and agrees well with CONTACT. (2) Under the elasto-plasticity effect, the maximum von Mises stresses and the maximum pressure are reduced and the shape of the contact patch changes from elliptic to egg shaped. Furthermore, the stick area becomes larger with both explicit and implicit methods and the maximum contact shear stress is slightly reduced. (3) The initial slope of the creep force characteristics obtained from the FE implicit solution under elastic conditions is slightly lower than CONTACT solution in the longitudinal and the lateral direction. Under elasto-plastic conditions, it is shown for a realistic example that the initial slopes of these curves are reduced by 15% in the longitudinal direction and by 13% in the lateral direction compared to the FE elastic solution. However, the relevance of these values should
be further investigated in more situations such as different wheel-rail couples in order to decide which reduction coefficients to use in vehicle dynamics software. (4) The FE implicit solution is sensitive to the energy dissipation related to the size of the time increment. Indeed, if the time increment is not small enough, oscillations at the rear of the contact area as well as inside the contact shear stress are observed which leads to stick/slip distribution widely different compared to CONTACT. Thus, the time increment must be carefully chosen by analyzing time increment sensitivity at the beginning of the study. Acknowledgments
The authors thank the Technological Research Institute RAILENIUM for their funding through the project CERVIFER and TATA STEEL for providing the material data References [1] F.W. Carter, On the action of a locomotive driving wheel, Proc. R. Soc. Lond. 112 (1926) 151-157. [2] K.L. Johnson, The effect of a tangential contact force upon the rolling motion of an elastic sphere on a plane, J. Appl. Mech., Trans ASME, 25 (1958) 338-340. [3] P.J. Vermeulen and K. L. Johnson, Contact of non-spherical bodies transmitting tangential forces, J. Appl. Mech., 31 (1964) 338-340. [4] K.L. Johnson, Contact Mechanics, Cambridge University Press, Cambridge, 1987. [5] J. J. Kalker, On the rolling contact of two elastic bodies in presence of dry friction, (Ph.D. thesis), Delft University of Technology, 1967. [6] J.J. Kalker, A fast algorithm for the simplified theory of rolling contact, Vehicle Syst. Dyn., 11 (1982) 1-13. [7] J.J. Kalker, Three-Dimensional Elastic Bodies in Rolling Contact, Kluwer Academic Publishers, Dordrecht, 1990. [8] C. Pritchard, Traction between rolling steel Surfaces - A Survey of Railway and Laboratory Experiments. Proceedings of the 7th Leeds/Lyon Symposium on Tribology, 1980. [9] A. Hobbs, A survey of creep, British Railways Research Dpt Report note DYN 52, 1967. [10] P.R. Nayak, S. Hariharan, R. Stern, R. Abilock and P.A. March, Friction and creep in rolling contact, Federal Railroad Administration, Dpt of transportation, 1970. [11] T.G. Pearce and K.A. Rose, Tangential force-creepage relationships in theory and practice, Contact mechanics and wear of rail/wheel systems, University of Waterloo Press, 1982. [12] E.A.H. Vollebregt, Numerical modeling of measured railway creep versus creep-force curves with CONTACT, Wear 314 (2014) 87-95. [13] A. Meierhofer, C. Hardwick, R. Lewis, K. Six and P. Dietmaier, Third body layer-experimental results and a model describing its influence on the traction coefficient, Wear 314 (2014) 148-154. [14] X. Zhao, Z. Li and R. Dollevoet, Solution of the wheel rail rolling contact in elasticity and elasto-plasticity using a transient finite element, International Conference on Contact Mechanics and Wear of Rail/Wheel Systems (CM2009), Firenze, ltaly, 2009. [15] X. Zhao and Z. Li, The solution of frictional wheel-rail rolling contact with a 3D transient finite element model: Validation and error analysis, Wear, 271 (2011) 444-452. [16] Z. Wen, L. Wu, W. Li, X. Jin and M. Zhu, Three-dimensional elastic–plastic stress analysis of wheel–rail rolling contact, Wear, 271 (2011) 426-436. [17] E.A.H. Vollebregt, User Guide for CONTACT, Vollebregt & Kalker's Rolling and Sliding Contact Model, Version13.1, VORtech, 2013. [18] Abaqus User Guide, Version 6.12, Dassault Systèmes Simulia Corp., 2012. [19] R. Courant, K. Friedrichs and H. Lewy, On the partial differential equations of mathematical physics, Math. Ann. 100 (1928) 32-74. [20] H.M. Hilber, T.J.R. Hughes and R.L. Taylor, Improved numerical dissipation for time integration algorithms in structural dynamics, Earthquake engineering and structural dynamics, 5 (1977) 283-292. [21] P. Aknin, H. Chollet, M. Sebes, Expro, Exploration de profils roue-rail, Version 4.0, INRETS, 2005. [22] P. Jehel, Modélisation numérique des phénomènes d'amortissement par dissipation d'énergie matérielle dans les structures de type portique en béton armé sous séisme (Ph.D. thesis), Ecole Normale Supérieure de Cachan. 2009. [23] M. Wiest, E. Kassac, W. Davesa, J.C.O. Nielsen and H. Ossberger, Assessment of methods for calculating contact pressure in wheel-rail/switch contact, Wear 265 (2008) 1439-1445.
Highlights
A 3D FE model is developed to study the frictional rolling wheel-rail contact
Explicit and implicit integration schemes are used and the solutions are compared
Implicit solution is sensitive to energy dissipation related to time increment size
The elasto-plasticity reduces the initial slopes of the creep force characteristics