Abstract
Over the past few years, a number of implicit/explicit finite element models have been introduced for the purpose of tackling the problems of wheel–rail interaction. Yet, most of those finite element models encounter common numerical difficulties. For instance, initial gaps/penetrations between two contact bodies, which easily occur when realistic wheel–rail profiles are accounted for, would trigger the problems of divergence in implicit finite element simulations. Also, redundant, insufficient or mismatched mesh refinements in the vicinity of areas in contact can lead to either prohibitive calculation expenses or inaccurate implicit/explicit finite element solutions. To address the abovementioned problems and to improve the performance of finite element simulations, a novel modelling strategy has been proposed. In this strategy, the three-dimensional explicit finite element analysis is seamlessly coupled with the two-dimensional geometrical contact analysis. The contact properties in the three-dimensional finite element analyses, such as the initial “Just-in-contact” point, the exact wheel local rolling radius, etc., which are usually a priori unknown, are determined using the two-dimensional geometrical contact model. As part of the coupling strategy, a technique has been developed for adaptive mesh refinement. The mesh and mesh density of wheel–rail finite element models change adaptively depending on the exact location of the contact areas and the local geometry of contact bodies. By this means, a good balance between the calculation efficiency and accuracy can be achieved. Last, but not least, the advantage of the coupling strategy has been demonstrated in studies on the relationship between the initial slips and the steady frictional rolling state. Finally, the results of the simulations are presented and discussed.
Keywords
Introduction
Rolling frictional contact between wheel and rail (W/R) is a highly non-linear problem involving large deformation, large rotation, material plasticity, contact, friction, etc. With the development of modern computing techniques and the availability of supercomputers, advances in the field of numerical simulations on such complex problems have been strongly boosted. Among the various numerical methods proposed,1–4 the finite element (FE) method is more widely used, by virtue of its striking versatility (i.e. accounting for arbitrary contact geometries, material plasticity, etc.). In general, based on the different features of solution procedures, the FE methods are classified into two main categorises,5,6 namely implicit and explicit.
Regarding the implicit FE method, a variety of models/tools have been created for different engineering purposes.7–12 For instance, Wiest et al. 7 performed implicit FE analyses to predict the normal pressure of W/R impact at a crossing nose. Telliskivi et al. 11 developed an implicit FE model to understand the complex behaviour of W/R interaction. The ratcheting performance of rail steels was evaluated by Pun et al. 8 The problem of normal contact was resolved using a quasi-static FE simulation, while the tangential shear stress distributions were calculated according to Carter’s theory. 4 The effect of wheel motions on the distribution of residual stresses was studied by Bijak-Żochowski and Marek. 10 Mandal and Dhanasekar 9 presented a novel sub-modelling approach to investigate the ratcheting failure of insulated rail joints. Based on the detailed stress/strain responses obtained from FE simulations, attempts were made by Ringsberg et al. 12 to predict the fatigue life of crack initiation on the rail surface. Ma et al. 13 introduced an implicit FE tool to qualitatively evaluate the performance of different rail pre-grinding strategies. However, due to the difficulties of convergence and the absence of dynamic effects,5,6 these implicit FE approaches were no longer able to meet the increasing expectations of FE-based contact models possessing higher degree of realism and better accuracy.
As an alternative problem-solving procedure (i.e. opposed to the implicit FE method), the explicit FE simulation proceeds without solving a large set of simultaneous equations at each time step and inverting the large matrix.5,6,14,15 This enables the explicit FE method to avoid certain difficulties of non-linear programming that the implicit method usually has. 16 Owing to such intrinsic advantages, the explicit FE approach is being increasingly welcomed for combatting the associated problems of W/R interaction. More recently, a series of representative three-dimensional (3D) explicit FE models17–20 have been presented. Zhao and Li 17 developed a 3D explicit FE model for W/R interaction. The model was successfully verified against CONTACT. An explicit FE tool was created to simulate the W/R contact–impact at rail insulated joint by Wen et al. 20 It was found that the variation of axle load had stronger effects on the impact event than other operational parameters. Pletz et al. 19 presented a dynamic wheel/crossing FE model to investigate the influence of operational parameters, such as axle loads, train speeds, material properties, etc., on the impact phenomena. The stress states and material responses under different levels of adhesions were analysed by Vo et al. 18 It was anticipated that the rail was highly prone to the damage resulting from the ratcheting fatigue of the material.
In summary, significant progress in the field of FE simulations on W/R interaction, from which a number of valuable insights are gained, has been made. Yet, there is still much work to be done. For example, the “gaps or penetrations” between wheel and rail interfaces cause the problems of divergence in the implicit analyses15,16 or even an unexpected failure in the explicit FE analyses. 21 Also, due to a priori unknown contact location, the general experience- or visualisation-based discretisation on the potential contact area could always lead to a redundant, insufficient or mismatched FE mesh. As a consequence, the efficiency and accuracy of FE simulations would be adversely affected.
To address these modelling challenges (explained later in detail) and to improve the performance of FE simulations, a new coupling strategy, that couples the two-dimensional geometrical (2D-Geo) contact analysis and the three-dimensional explicit finite element analysis, has been proposed. The idea of this strategy is inspired by the approach used in CONTACT (i.e. a well-established computational programme developed by Professor Kalker 22 and powered by VORtech Computing 3 ). Prior to the 3D-FE analysis, the 2D-Geo contact analysis is functioning effectively as an adaptive guidance for the identification and discretisation of the potential area in contact. The challenge of the coupling strategy lies in the programming efforts on how to build up the 2D-Geo and 3D-FE models as well as how to make the interface of two models work effectively. Based on the simulation results, it has been noticed that the computational problems, such as gaps/penetrations, mesh refinement, unexpected initial slips, etc., have been resolved successfully. The goal of improving the performance of FE simulations has thus been achieved with the aid of this coupling strategy.
The outline of this paper is as follows. In the next section, full attention is focused on the general descriptions of the 3D-FE model. Also, the details of FE modelling challenges, that prohibit the analysts from attaining accurate contact solutions, are presented. The strategy (referred to as ‘enhanced explicit FE-based coupling strategy’, abbreviated as ‘eFE-CS’, hereinafter), which couples the 2D-Geo model together with that of 3D-FE, is described in the Coupling strategy section. Following that, the effectiveness and advantages of ‘eFE-CS’ strategy are demonstrated in the FE results and discussion section. Finally, concluding remarks are drawn.
W/R 3D-FE model
In this section, the FE model for the analysis of W/R interaction is presented. The two counterparts investigated here are the standard S1002 wheel 23 (EN13715-S1002/h28/e32.5/6.7%) with a nominal rolling radius of 460 mm and the standard UIC 54E1 rail. Here, ‘h28’ refers to the flange height of 28 mm; ‘e32.5’ means the flange thickness of 32.5 mm; ‘6.7%’ is the reverse slope. The drawing of wheel cross section is adopted from International Union of Railways. 24 The inner gauge of the wheel-set is 1360 mm and the axle length is 2200 mm. The track gauge is 1435 mm. Also, the cant angle of 1/40 is considered in the model. Note that the model can be easily adjusted to account for other different W/R profiles.
Discretised FE model
A half W/R FE model shown in Figure 1(a) and (b), where the rail is modelled with restriction to an overall length of 1.8 m, is adopted by taking advantage of the symmetrical characteristic of the track and the wheel-set. Such a half W/R FE model is similar with and inspired by the ones described in the literature.19,20 The wheel is set to roll from the origin of the global coordinate system FE model of the W/R interaction: (a) schematic graph and (b) front view.
Only the regions where the wheel travels are discretised with dense mesh, leaving the remaining regions with coarse mesh (see Figure 1(a) and (b)). A solution area is introduced and positioned in the middle of the dense meshed area. Here, the solution area is defined as a region to extract and analyse the contact properties, such as the contact patch, normal pressure, shear stress, etc. In this region, the mesh size (1 mm) is approximately two times smaller than the dense meshed area (2 mm) for the purpose of capturing the associated high stress/strain gradients.
Material properties and mechanical parameters.
Note: “a”: Linear elastic material.
To better simulate the process of W/R dynamic contact, an implicit-to-explicit sequential solving procedure 21 is adopted. In this procedure, the implicit solver (ANSYS Mechanical 15 ) and explicit solver (ANSYS LS-DYNA 21 ) work in pairs. First, the equilibrium state of the preloaded structure (e.g., under the axle-load of 100 kN) is determined with ANSYS Mechanical. The displacement results of the implicit analysis are used to do a stress initialisation for the subsequent transient analysis. Then, the processes of rolling frictional contact begin at time zero with a stable preloaded structure. 15 The responses of dynamic contact are further simulated with ANSYS LS-DYNA following the scheme of central difference time integration. 21
For such a typical FE analysis of dynamic contact, the basic process consists of three tasks: (1) Build the model, including prescribing the initial location of W/R, defining correct boundary conditions, preforming mesh refinement, etc. (2) Apply the loads and run the simulations, involving traction application, contact definition and determination of calculation time step size; (3) Review the results, referring to the visualisation of contact properties, such as surface normal pressure, shear stresses within the actual contact patches, subsurface stress/strain responses, etc.
Challenges of FE analysis
As the contact always occurs at a priori unknown area, a series of obstacles will thus be encountered when performing FE-based contact simulations. In this section, the details of those obstacles are presented and discussed.
Geometrical gaps/penetrations
Figure 2 shows the ‘gaps/penetrations’ between the contacting bodies of W/R. Both the gaps and penetrations are undesired. For example, in the implicit FE analyses, too large gap may lead to the wrong identification of contact pairs, while prominent penetrations may cause the overestimation of contact forces. Consequently, it causes an unexpected failure (usually the error of divergence) of the FE simulation. For explicit FE analyses, no penetrations are allowed and the gaps can increase the calculation cost significantly for eliminating the effect of initial disturbances.
Schematic of initial contact status: (a) penetration, (b) gap and (c) just-in-contact.
Although there are several options (i.e. tuning the real constants and key options of contact elements independently or in combination) suggested in ANSYS to adjust the initial contact conditions, they do not work well when the values of gaps/penetrations are large (e.g. in the order of millimetres or even larger). Also, it is quite challenging to manually bring the two W/R contacting bodies into the positions of “Just-in-contact” as shown in Figure 2(c). Here the term of “Just-in-contact” refers to a contact positioning between wheel and rail, where the two bodies touch each other without or with a tolerable contact gap or penetration. Moreover, due to the complexity of W/R realistic contact geometry, it is very difficult to determine where the first point of contact will occur. Thus, the gaps/penetrations appear to be inevitable and troubling.
Preload application
As shown in Figure 3, the wheel is prescribed to travel over the rail from the initial condition of pure rolling state to the steady state of partial slip. Here, the steady state
22
means that the slip-adhesion phenomenon is independent of time. To achieve the initial state of pure rolling without slip, the velocity Schematic of the wheel kinematics (note: “C” denotes the wheel centre-of-mass).

As the “Just-in-contact” point is a priori unknown, the real local wheel rolling radius
In addition, to reproduce the local phenomenon of partial slip numerically, the appropriate traction torque
Mesh refinement
Given that the refined element size inside the contact patch should be as small as 1000 times over the dimension of W/R components, an economic, adaptive and reliable mesh generation has always been a challenge for the analysts.
However, there are no clear rules proposed for determining the size of the refined potential contact area up to now. The most commonly adopted approach is the trial-and-error method mainly depending on the researchers’ experience or visualisation. Using this method, the refined potential contact area can be easily overestimated or underestimated (see Figure 4(a) and (b)). In particular cases, the refined potential contact areas might be mismatched or deviated from each other (see Figure 4(c)). Thus, the trial-and-error method is highly prone to the inaccurate or undesired contact solutions of W/R interaction.
Schematic of the mesh refinement at the potential contact area: (a) Redundant, (b) Insufficient, (c) Mismatched and (d) Idea.
In addition, if the relative contact locations between the wheel and the rail vary along the track, the corresponding mesh refinement is in demand to be altered spontaneously. This brings about an even higher requirement on the flexibility of the mesh refinement approach.
Coupling strategy
To deal with the aforementioned challenges of FE analyses, a novel coupling strategy, that combines 2D-Geo contact analysis and 3D-FE analysis, is developed (see Figure 5). The purpose of the 2D-Geo simulation is to detect and determine the initial “Just-in-contact” location (CP), the contact clearance and the corresponding local wheel rolling radius at the point of contact. The obtained contact information is used as the adaptive guidance for the 3D-FE analysis (see W/R 3D-FE model section).
Working mechanism of the coupling strategy.
To implement this coupling strategy, the data exchange between 2D-Geo and 3D-FE models is performed by a small external routine written in the MATLAB environment. Once the processes of 3D-FE modelling and preloading are accomplished, the dynamic simulation is performed to study the behaviour of W/R contact. The details of the 2D-Geo analysis and the interfacing scheme of the two models are briefly described in this section.
2D-Geo analysis
In order to better illustrate the 2D-Geo contact model, a rigid wheel-set is positioned over a rigid track shown in Figure 6(a). The global coordinate system (a) W/R 3D coordinate systems and (b) “Just-in-contact” equilibrium condition of the wheel–rail at a lateral displacement of −10 mm. (a) “Just-in-contact” equilibrium location derived from the 2D-Geo contact analysis and (b) the cross sectional view of the W/R dynamic contact FE model.

The current 2D-Geo contact model is further developed on the basis of the previous work. 13 The initial contact points, where the two particles on the un-deformed wheel and rail coincide with each other, are determined under a given lateral displacement of the wheel-set. The solution procedure of 2D-Geo model is non-iterative by taking advantage of efficient matrix operations in MATLAB, which means that no inner and/or outer loop iterations (e.g. ‘for’ or ‘while’ loops) of exploring the “Just-in-contact” equilibrium conditions are performed. By this means, a significant increase in the calculation efficiency (around 10 s) is accomplished, as opposed to the conventional iterative fashion (e.g. 2–3 min in Ma et al. 13 ). More information about the contact searching schemes is given in the literature.13,27,28 Figure 6(b) shows a typical example of the 2D-Geo contact analysis, in which the wheel-set is located at the “Just-in-contact” equilibrium position with a lateral displacement of −10 mm.
Coupled interface
Using the 2D-Geo model described above, the interface and outcome of the ‘eFE-CS’ strategy are presented below.
“Zero” gaps/penetrations and obtained Rw
Figure 7(a) shows the results of the 2D-Geo contact analysis, through which the 2D wheel-set is positioned on the “Just-in-contact” point of the track. Depending on the lateral displacement and roll angle of the wheel-set determined by the 2D-Geo analysis, the wheel position in the 3D W/R FE model (see Figure 7(b)) has been properly adjusted. The resulting gaps or penetrations between wheel and rail can be reduced to the order of a micrometre or even less. This ensures a successful FE simulation of W/R interaction by taking advantage of the suggested options in ANSYS (see Geometrical gaps/penetrations section).
Besides, the actual rolling radius
New adaptive mesh refinement
Considering that the refined potential contact area should completely encompass the actual contact area3,22 (see Figure 8(b)), the calculated contact clearance from 2D-Geo analysis is used as the adaptive guideline for the mesh refinement of W/R interface. Here, the term of “Contact clearance” Schematic representation of the new adaptive refining approach: (a) global view of contact clearance; (b) schematic graph of the potential and actual contact areas; (c) close-up view on the variation of contact clearance at the vicinity of initial contact points and (d) close-up view of the refined potential contact region.
For the sake of identifying the best contact clearance, it is expected to follow a guideline as follows:
The contact clearance The refined potential contact area determined by the contact clearance
Based on the selected value of the contact clearance, the mesh refinement process in these regions will be initiated. Using the proposed Nested Transition Mapped mesh refining approach,
29
the elements in the solution area are able to be refined as small as 1.0 mm × 1.0 mm × 1.0 mm (as shown in Figures 8(d) and 10), while for the remaining out-of-contact region (coarse mesh area), less attentions will be paid.
Flexibility of the new mesh refining approach with respect to. different lateral displacements: (a) 5.5 mm; (b) 2.0 mm; (c) −3.5 mm and (d) −5.5 mm. The resulting contact pair: (a) isometric view and (b) side view.

Flexibility of the new mesh refining approach
To automate the common tasks and implement the coupling strategy more efficiently, both the 2D-Geo and 3D-FE models have been coupled and parametrised using the scripting language of MATLAB and Parametric Design Language of ANSYS (APDL). The scripted and later packaged ‘eFE-CS’ modelling strategy offers great convenience and flexibility for the day-to-day analyses.
One added important feature of this adaptive mesh refining method is shown in Figure 9, where the patterns of these local mesh refinements change flexibly with respect to different lateral displacements of the wheel-set (e.g. ranging from −5.5 mm to 5.5 mm). For the cases of non-zero lateral displacements, the rails remain symmetric in relation to the track centre, while only the wheel and contact conditions vary. The variations of geometrical contact properties associated (i.e. locations and dimensions of refined potential contact areas) are predicted by a priori 2D-Geo contact analysis (see Figure 6(b)). Once the FE models (i.e. with half W/R considered) are created, they are used to study the dynamic behaviour of W/R interaction. Here, the applicability of this half W/R FE modelling approach (as shown in Figure 1) to the cases of different lateral displacements are explained as:
A complete wheel-set and two rails are considered in the 2D-Geo contact model (see Figure 6). The roll angles of the wheel resulting from the variation of lateral displacements have been explicitly taken into account. A range of lateral displacements (i.e. (−5.5 −5.5) mm) cover most sections of railway track (e.g. tangent, large radius curves, etc.),
30
where the yaw angles (also called attack angle) are usually small (less than 2°
31
). Moreover, being aware that the wheel travels over a relatively short distance of 0.52 m, which is rather short in comparison to the wavelength of yaw/Klingel motion (i.e. usually in the order of 10–100 m), the variation of yaw angle would be relatively small. It is thus assumed that the neglect of yaw motion in the model presented might be acceptable. A comprehensive model verification against CONTACT
3
has been performed by the authors recently.
32
For the cases of different lateral displacements (i.e. (−5.5 −5.5) mm), the model verification shows that the half W/R FE model can produce as good results as CONTACT.
In summary, the half W/R FE model enhanced with 2D-Geo contact analysis is flexible enough to be used in the cases of different lateral displacements. However, the further applicability of half W/R FE model to the sharp curves (i.e. curve radius smaller than 500 m 33 ), in which the yaw angle is indispensable and plays an important role, will be explained later in the Discussion: Pros and cons of ‘eFE-CS’ model section.
Resulting contact pairs
Figure 10(a) and (b) shows a contact pair, consisting of two master and slave segments for explicit FE analysis (TARGET 170 and CONTACT 173 elements for implicit FE analysis). It can be seen that the size of the potential contact area is controlled by the magnitude of the refined width Δ
It has been reported in the literature34,21 that contact analysis will add a significant computational cost for the overall solution, even when the ratio between the number of contact points and the number of elements is small. To limit the level of calculation expense as low as possible, the smallest and localised contacting zones are always desired to maintain the most efficient solution. Also, as discussed before, the refined potential contact areas should be large enough to encompass all the necessary contacts to maintain the accuracy of the solution. According to the coupling strategy, the refined width Δ
FE results and discussion
To demonstrate the effectiveness of the ‘eFE-CS’ strategy, a series of FE simulations have been performed. The influence of the size of potential contact areas and the local wheel rolling radius on the performance of ‘eFE-CS’ strategy are to be analysed. As the wrongly estimated rolling radius
Contact clearance
Figure 11 shows the effect of contact clearances ( Variation of refined widths with respect to different contact clearances: (a) relations between contact clearance Comparison of FE models with respect to three different refined widths. Number of elements in FE model. Calculation expense.
Since the actual contact areas would not vary much (i.e. the same axle load is applied in both implicit and explicit FE analyses), the effect of varying contact clearances
Figure 12 shows the variation of normal pressure corresponding to different contact clearances. A noticeable discrepancy of the maximum contact pressure is observed between Cases I (1216 MPa) and II (1041 MPa). Such a discrepancy is explainable, since the mesh sizes of FE models at the origin of wheel rotation (not the solution area, see Figure 1) are relatively large (around 2 mm). Also, considering that W/R contact occurs at a highly stressed area (i.e. as small as 150 mm2, see Ma et al.
35
), these large mesh sizes are hardly possible to capture the stress/strain gradients accurately. It can thus readily trigger the tolerable discrepancies of the maximum contact pressure. (Notation: the contact responses extracted from the solution area, where the mesh size is as small as 1 mm, are preferred.)
Variation of actual contact patches and contact pressure with respect to different refined widths: (a) Case I: 
With respect to Case III (
Moreover, it can be observed from Table 2 that the amount of the discretised elements as well as the calculation expense reduce significantly with the decrease of the contact clearance
To assess the applicability of the selected contact clearance (0.17 mm) to varying axle loads, two more case studies with axle load of 120 kN and 140 kN have been performed. Figure 13 shows the comparison of actual contact patches with respect to different axle loads. It can be seen that there is a steady increase in both the maximum contact pressure and the actual contact area in relation to the increasing axle loads. The width of refined potential contact areas determined by the contact clearance of 0.17 mm is robust enough to encompass all the real contact areas, which are resulting from the applied axle loads up to 40% (140 kN) higher than the normal one (100 kN).
Comparison of actual contact patches and contact pressure at different axle loads: (a) 100 kN. (b) 120 kN and (c) 140 kN.
Similarly, the applicability of the selected contact clearance of 0.17 mm to the varying lateral shifts of the wheel-set (see Figure 9) has been evaluated. The results are shown in Figure 14(a) to (c). It can be seen that although the contact clearance of 0.17 mm can satisfy the case with the lateral shift of 5.5 mm (Figure 14(c)), it is insufficient for the case of −5.5 mm (Figure 14(a)), where the real contact patch tends to fall out of the refined potential contact area. Given a larger contact clearance Comparison of actual contact patches and contact pressure at different lateral shifts: (a) −5.5 mm (
In summary, a contact clearance of 0.17 mm performs well at zero lateral displacement of the wheel-set. It is robust enough to apply in the cases of higher axle loads. One can further affirm that such a contact clearance is applicable to the case of different operational conditions (i.e. train velocities, friction coefficients, traction coefficients, etc.), since their influences on the size of real contact patches are insignificant in comparison to the increased axle loads. The applicability of contact clearance 0.17 mm to these varying operational patterns has been verified in the authors’ recent work.32,35
However, with respect to the cases of large (e.g.
In extreme cases (e.g. severe worn contact geometries), if the suggested contact clearance (either 0.17 mm or 0.35 mm) cannot guarantee good results, it is recommended to follow the above presented parametric study procedure to make a good decision. Normally within two (maximally three) times of targeted parametric examinations of contact clearances, the best refined width of the potential contact area can be found.
Travelling distance
Apart from the contact clearance
Figure 15(a) shows the results of contact forces obtained from implicit-to-explicit FE analyses. It can be observed that the W/R mutual reactions are equal in magnitude and opposite in direction. Such observations are complementary to the classical contact theory as demonstrated in McMichael.
36
(a) Responses of W/R contact forces with a travelling distance (
The solid lines represent the contact response in the dense meshed area, while the diamond dotted lines denote that inside the solution area. It is observed that the vertical contact forces
In addition, the longitudinal frictional forces of
In comparison with the dense meshed region (see Figure 15(a)), all the force components exhibit a “sudden perturbation”, as the wheel rolls and approaches the vicinity of the solution area (80 mm). The causes of these “sudden perturbations” have been studied extensively and found (see supplement material)
35
to be the difference of the associated mesh sizes, which vary from 2 mm (the dense meshed area) to 1 mm (the solution area). For this reason, the magnitude of nodal contact stiffness
Accordingly, the change of contact stiffness manifests itself in the variation of contact forces
Here,
To sum up, the variation of mesh sizes implies the “sudden perturbations” of contact forces as shown in Figure 15(a). Also, it has been reported in Ma et al.
35
that the problems of “sudden perturbations” can be addressed using an optimal penalty scale factor
Comparison of FE models with respect to different travelling distances.
Implicit-to-explicit sequential solution.
Number of elements.
Calculation expense.
Figure 16(a) shows three of the five representative refined FE models with the travelling distance of 240, 520 and 800 mm, respectively. A close-up view on the performance of these varying travelling distances is taken by a supplemental and independent check on the resulting contact forces (as suggested by McMichael
36
). The results shown in Figure 16(b) indicate that all the resulting vertical forces are varying around the axle load (100 kN) exerted on the wheel axle. It is clear that the oscillations of contact forces (“sudden perturbations” as depicted with diamond dotted lines) occur for all the five cases of different travelling distances.
Influence of varying travelling distances on the W/R dynamic contact responses: (a) representative FE models of the different travelling distances and (b) vertical contact force with respect to different travelling distances.
Similarly, the “saw-toothed” force oscillations, which have been explained previously, are observed in Figure 16(b). Also, it can be seen that all the FE models take around 150 mm to completely get rid of the initial contact disturbances under the specified operational conditions as shown in Table 1.
Figure 17 shows the normal contact pressure that is extracted at the instant when the wheel travels over the middle of the solution area. It can be seen that both the magnitude and the distribution of the contact pressure vary slightly with respect to different travelling distances.
Effect of varying travelling distances on contact pressure.
In view of the increasing calculation expense (as listed in Table 3) as well as the varying contact forces and pressures, it is practically necessary to build up a guideline/criterion to select the best travelling distance. Aiming to maintain the good compromise between the calculation efficiency and accuracy, the guideline is formulated as a recommendation:
The overall travelling distance should be kept shorter enough to save the calculation expenses. A distance of 150 mm before the solution area should be reserved for damping out the “saw-toothed oscillations”.
Following the aforementioned guideline, the most suitable travelling distance of 520 mm is suggested. This suggested that the travelling distance
Initial slip
To study the influence of initial slips
Figure 18 shows the variation of longitudinal frictional forces with respect to different initial slips Variation of frictional forces with respect to different initial slips. Effect of varying the initial slips on the surface shear stress.

With the increase of initial slips
Similar patterns of shear stress distribution occur in the other two cases, Cases IV (
When the value of
Such phenomena could be explained by the wheel kinematics as shown in Figure 20, which depicts the relationship between the frictional forces and the slip velocities. It can be seen that the directions of frictional forces are in opposite to the ones of slip velocities. Thus, the magnitude of frictional forces changes from negative to positive in the Cases V and VI, where the initial slips are positive as opposed to the others.
Schematic of the wheel kinematics: (a) negative initial slip 
Figure 20(c) schematically shows the traction curve.
25
It is observed that when the magnitude of the slip (also interpreted as creep-age) exceeds certain level, the regime of full slip is entered. The frictional force thus reaches its saturation value (i.e. traction bound of 50 kN). Here, the traction bound, which amounts to the product of frictional coefficient
For the case of initial slip
To summarise, using the proposed ‘eFE-CS’ approach, the local wheel rolling radius
Discussion: Pros and cons of ‘eFE-CS’ model
From the results of the FE simulations presented, it can be noticed that the proposed ‘eFE-CS’ strategy is promising enough for improving the performance of FE simulations on W/R interaction. The advantages it holds can be categorised into four groups:
All the modelling challenges (see Challenges of FE analysis section), from which most of 3D-FE models often suffer, can be readily addressed. Also, a good compromise between the calculation accuracy and efficiency is maintained. More detailed normal and tangential contact solutions than those of 2D W/R contact analysis37,38 are generated. These stress/strain responses are expected to contribute positively on other advanced applications (e.g. prediction of wear/rolling contact fatigue, profile design/optimisation, etc.). The dynamic effects, which are often neglected in the static or quasi-static 3D FE models
13
(contact constrains are enforced with Lagrange Multiplier method, Arbitrary Lagrangian Eulerian method, etc.), are taken into account in the present model (Penalty method). A high degree of realism is attained. The strategy shows good applicability to the problem of W/R interaction. As there are no special restrictions imposed upon its use, it is recommended to use for addressing other contact/impact problems (e.g. gear, bearing, metal forming, etc.) that have complex contact geometries.
Although the developed ‘eFE-CS’ model has many advantages, there is room remained to further increase its degree of realism and accuracy. For instance,
To enhance the steering capability: In the present ‘eFE-CS’ model, the steering capability of the wheel is limited, since the yaw motion is not allowed and only roll motion is considered. For this reason, it would be difficult to study the cases of W/R interaction at sharp curves (curve radius of 500 m or even small). In the future work, it is motivated to integrate several torsional control units onto the wheel for the purpose of attaining a sufficient steering capability. To account for the accumulation of wheel passages: It is known that the degradation of W/R interface is due to cyclic contact loadings. To better predict the degradation process, the modelling of multiple wheel passages
39
is increasingly demanded. To study the effect of second wheel of the axle: Using a half W/R FE model, the second wheel of the axle will have a mirrored roll angle and contact conditions, which correspond less well to the reality. As a consequence, the accuracy of the FE simulations would be adversely affected. For a more quantitative explanation in terms of how the mirrored contact conditions affect the contact solutions, rigorous verifications on the half W/R FE model against the complete wheel-set/track FE model are needed to be performed in the future work.
Conclusions
Aiming to improve the performance of FE analysis on the wheel–rail interaction, a new ‘eFE-CS’ strategy, that couples the 3D-FE model and the 2D-Geo contact model, has been presented. Prior to the 3D-FE simulations, the 2D-Geo contact analysis is performed to better define the contact properties that are in demand for FE analyses. Both the advantages and disadvantages of the developed ‘eFE-CS’ model have been discussed. Based on the simulation results and discussions, the following conclusions are drawn:
Following the coupling procedure, the convergence problems or even a failure of the FE simulations in the presence of gaps or penetrations have been addressed by correctly placing the wheel onto a position of “Just-in-contact” over the rail. In order to avoid the issues of the redundant, insufficient or mismatched mesh refinement in the vicinity of the contact region, which will lead to either the prohibitive calculation expense or the inaccurate solution, an adaptive mesh refining technique based on the contact clearance has been proposed. Using this technique, the solution of W/R rolling frictional contact is able to be maintained with good accuracy and efficiency. The flexibility of this coupling approach has also been demonstrated by solving the cases of different contact locations due to various lateral shifts of the wheel-set. An unexpected initial slip can be introduced by mismatching the angular and translational velocity of the wheel, owing to the wrongly estimated local wheel rolling radius. Such an initial slip can cost unnecessary calculation effort to achieve the steady frictional rolling state or even lead to wrong tangential solution of W/R interaction in extreme cases. Using the coupling strategy, the value of undesired initial slip can be easily minimised. Accordingly, the accurate and efficient FE solution can be guaranteed.
In conclusion, the proposed coupling strategy allows the proper prescription of the contact properties, which can enhance the performance of the FE analyses of W/R interaction significantly. Such a strategy is promising enough to promote the application of the FE approaches on the contact problems with more complex and irregular contact geometries (e.g. the mild or severe worn W/R interfaces, metal forming, gear, bearings, etc.).
Footnotes
Acknowledgements
The comments and suggestions from Prof. Rolf Dollevoet on the manuscript are gratefully acknowledged. Special thanks go to Dr Hongxia Zhou for critically reading this manuscript and giving helpful suggestions. The authors are also very grateful to all the reviewers for their thorough reading of the manuscript and for their constructive comments, which have helped to improve the manuscript.
Declaration of Conflicting Interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Yuewei Ma thanks the CSC (China Scholarship Council) for their financial support.
Supplementary Material
Supplementary material is available for this article online.
