Wheel–rail impact loads and axle bending stress simulated for generic distributions and shapes of discrete wheel tread damage

Wheel–rail impact loads generated by discrete wheel tread irregularities may result in high dynamic bending stresses in the wheelset axle, leading to a decrease in component life and an elevated risk for fatigue failure. In this paper, a versatile and cost-eﬃcient method to simulate the vertical dynamic interaction between a wheelset and railway track, accounting for generic distributions and shapes of wheel tread damage, is presented. The wheelset (comprising two wheels, axle and any attached equipment for braking and power transmission) and track with two discretely supported rails are described by three-dimensional ﬁnite element (FE) models. The coupling between the two wheel–rail contacts (one on each wheel) via the wheelset axle and via the sleepers is considered. The simulation of dynamic vehicle–track interaction is carried out in the time domain using a convolution integral approach, while the non-linear wheel–rail normal contact is solved using Kalker’s variational method. Wheelset designs that are non-symmetric with respect to the centre of the axle, track support conditions that are non-symmetric with respect to the centre of the track, as well as non-symmetric distributions of tread damage on the two wheels (or irregularities on the two rails) can be studied. Time-variant stresses are computed for the locations in the wheelset axle which are prone to fatigue. Based on Green’s functions for stress established using the wheelset FE model, this is achieved in a post-processing step. An extensive parametric study has been performed where wheel–rail impact loads and axle stresses have been computed for different distributions and sizes of tread damage as well as for different train speeds.


Introduction
Wheel tread damage leading to high magnitudes of vertical dynamic wheel-rail contact forces is a major cause of train delays in the Swedish railway network, in particular during the coldest months of the year. According to regulations, vehicles generating contact forces exceeding the limit value for allowed wheel-rail impact loads must be taken out of service. This may lead to severe disruptions in the schedule of railway operations and unexpected costs for the train operators. Wheelrail impact loads induced by wheel tread irregularities (out-of-roundness), such as wheel flats and rolling contact fatigue (RCF) clusters, result in elevated stress levels in wheels, axles and bearings, thus increasing the risk of fatigue failures [1] . In addition, increased loads due to wheel tread defects may shorten the life of track components, such as rails and sleepers, resulting in higher costs for maintenance of infrastructure. Wheel tread irregularities also lead to increased levels of rolling noise, impact noise and ground-borne vibration [2 , 3] .
As a consequence, careful planning of wheelset maintenance is of high importance for the train operators. Traditionally, regular time intervals for the inspection of wheelsets are scheduled depending on the mileage covered by the vehicle, see e.g. Ref. [4] . However, a more flexible maintenance planning could reduce overall maintenance costs, avoid non-necessary replacements of components and increase the availability of the fleet. Indeed, many operators are currently investing in procedures for condition-based maintenance [5] . By acquiring better knowledge of asset degradation and limits for reliable utilisation, this may be achieved without reducing the safety level in operations.
Condition-based maintenance is referred to as a proactive strategy [6] , where the maintenance intervals are optimised according to the observed conditions of the asset. As an example, the status of a given wheel tread surface can be monitored regularly by measuring the dynamic loads that are generated while the wheelset is passing over wheel impact load detectors distributed in the railway network [7] . Information on the status of a wheelset, and of the track, can also be gathered by monitoring bending and torsion stresses at critical sections of the axle [5] . In this way, maintenance can be postponed as long as the monitored parameters are within the "safe" range, which is determined by regulations or by assessment according to company practises [4] .
Condition-based wheelset maintenance can be enhanced by improving the knowledge of the influence of different forms of running surface damage on wheel and rail on component fatigue life. This paper will present a computational procedure to predict dynamic wheel-rail contact forces and wheelset stresses generated by generic forms (shapes, sizes and distributions) of wheel/rail irregularities. Rail and wheel acoustic roughness, wheel out-of-roundness (polygonalisation), rail squats, discrete wheel tread and rail head defects, as well as RCF clusters, are amongst the possible forms of irregularities that may be studied. In this paper, a mathematical description of a breakout of material from the wheel tread will be implemented to demonstrate the approach. A broad range of operational conditions can be studied, facilitating investigations of the influence of train speed, axle load, track and wheelset design, track stiffness irregularities, primary suspension characteristics, etc.
Several models and algorithms for the computation of wheel-rail contact forces have been presented in the literature, see e.g. the surveys in [3 , 8 , 9] . A general approach for solving the rolling contact problem, where discrete irregularities on wheels and rails are considered, is to use the FE method in the time domain. However, FE computations in the time domain are time demanding, in particular if a dense mesh of the wheel is employed and a long stretch of track needs to be considered. Computational times can be substantially reduced if the analysis is performed using a frequency-domain model. However, such models need to be completely linear. In studies where the variation in contact force is large with respect to the static wheel load, such as for the impact loads generated by the tread damage studied in the present work, a non-linear contact model is necessary, which requires calculations in the time domain [3] .
A common strategy to limit computational cost in the time domain is to introduce modelling simplifications. Most models for simulation of vertical dynamic interaction between vehicle and track assume that the vehicle and track properties, and also the excitation, are symmetric with respect to the track centre, see e.g. Refs. [10][11][12] . In this way, it is possible to reduce the size of the model while only accounting for the dynamics of a set of wheels (on one side of the train) interacting with one rail. An alternative, attractive approach to reduce the computational cost of a time-domain model is to model the wheelset and track by means of their impulse response functions, also known as Green's functions. These are pre-calculated before carrying out the simulation of dynamic vehicle/track interaction. This approach can be traced back to Heckl's proposal for a railway simulation program [13] .
In the current work, the model from Ref. [10] using such a Green's functions approach has been extended by accounting for both wheels of the wheelset, both rails of the track and the dynamic coupling between the two wheel-rail contacts that takes place via the wheelset axle and via the sleepers. This means that a non-symmetric wheelset design with for example a laterally offset gearbox assembly can be considered. Further, generic combinations of non-symmetric tread damage on the two wheel treads in the same wheelset can be studied, as well as non-symmetric track support conditions. Despite the extension of the model, the computational cost of the time-domain calculation is low due to the usage of pre-calculated Green's functions for the dynamics of the wheelset and the track. The adopted computational procedure uses an in-house implementation of Kalker's variational method [14] to solve the vertical wheel-rail contact forces when accounting for a generic geometry (shape, size and distribution) of wheel tread and/or rail head defects that are superimposed on the nominal contact surfaces of the wheel and rail. For each wheel-rail contact, the solution of the contact problem provides the time history of the vertical contact displacements of wheel and rail, as well as the time history of the vertical contact force. The calculated time variations of contact forces are then used to compute the dynamic stresses generated in the wheelset.
Traditionally, a wheelset is designed by the use of specific standards, such as EN 13,104 [15] for powered axles and EN 13,979-1 [16] for wheels. However, these standards are based on static calculation procedures where magnification factors are used to account for the dynamic variations in loads. EN 13,104 can be used to calculate fatigue stresses at different  sections of the axle (based on axle design, information on vehicle dimensions, operating conditions, and braking or traction  system) and comparisons can be made with allowed stress levels for specific axle steel grades at specified axle zones. Similarly, EN 13,979-1 details how design loads based on the static axle load can be used to calculate pertinent fatigue stresses for wheels (by use of the FE method) that can be compared to allowed stress levels. Common for these standards is that the implemented design forces and moments acting on the wheelset are independent of the actual operating conditions and do not account for the degradation of the wheelset during operation.
By applying the computational procedure proposed in this paper, the calculated vertical contact force variations resulting from prescribed surface irregularities of wheel tread can be used to compute the stress variations at critical positions in a wheelset. A similar approach was adopted in Refs. [17 , 18] , where the dynamic stresses were obtained by a processing of the displacements associated with the wheelset flexibility and rigid body modes. The displacements were calculated based on the mode shapes of the wheelset. However, only axle points sufficiently far away from stress concentrations could be studied in [17] .
Calculated stress variations can be compared to those determined by the vertical loads given in the standards. Based on such comparisons, or fatigue life calculations using a suitable high-cycle fatigue criterion or a crack growth algorithm, the planned maintenance can be scheduled to an earlier or later occasion. Priority can be given to trains where the monitored stresses of some components are more critical. This allows for a more efficient use of resources for the train operator, as well as more flexibility in the usage of the fleet.
In the following sections of the paper, the applied theory for the simulation of vertical dynamic wheelset-track interaction, including the implementation of the cross-coupling between the two wheel-rail contact points and the stress calculation procedure, is presented. Then the wheelset and the track models are described. Finally, some examples of calculated contact forces and axle stresses for different running conditions and distributions of tread damage are discussed to demonstrate the approach.

Theory
The present section describes the approach that has been implemented to simulate the vertical dynamic interaction between a wheelset and railway track. The models of wheelset and track, as well as the input data that have been adopted, will be described in more detail in Chapters 3 and 4. Most of the calculations described here are performed using the in-house MATLAB software WERAN for simulation of dynamic vehicle-track interaction based on the Green's functions approach, originally presented in Refs. [10 , 11] and later extended in Ref. [19] to also include lateral wheel/rail interaction. A validation of the modelling approach against another model for simulation of dynamic vehicle-track interaction (DIFF) was presented in Refs. [20 , 21] . DIFF in turn has been validated against field measurements [22] .
The novel contribution of this paper is the extension of the model in Ref. [10] to also make possible simulations of wheelset designs that are non-symmetric with respect to the centre of the axle, non-symmetric distributions of tread damage on the two wheels (or irregularities on the two rails), as well as track support conditions that are non-symmetric with respect to the centre of the track. The model accounts for the coupling between the two wheel-rail contacts (one on each wheel) via the wheelset axle and via the sleepers. This means that two contact problems are simulated simultaneously, one for each nominal wheel-rail contact point on the same wheelset. A consistent approach to calculate stresses in the wheelset based on a convolution integral formulation is also presented. Some examples will be demonstrated in Chapter 5.

Green's functions for wheelset and track
In this paper, Green ś functions are established from frequency response functions determined by three-dimensional (3D) FE models of wheelset and track developed in the commercial FE software ABAQUS 2018 [23] . The complex-valued (including information about both magnitude and phase) frequency response functions are calculated using a subspace-based steadystate dynamic analysis [23] . For the wheelset, a single vertical harmonic load is applied on the nominal rolling contact point on one of the two wheels and the resulting displacement responses are extracted at both nominal contact points to provide the direct and cross receptances of the wheelset at uniform frequency intervals. For a non-symmetric wheelset design, the direct receptances are different for the two wheels (meaning that the procedure needs to be duplicated for both wheels), while the cross-receptances are the same due to reciprocity.
Once the receptance R W i, j ( f ) has been calculated at point i of the wheelset for a load applied at point j , the corresponding Green's function G W i, j (t) is calculated by means of an inverse Fourier transform symbolically expressed as: Since the rotation of the wheelset is not accounted for here, the Green's function G W i,i is equal to the impulse response of the radial displacement at point i due to radial excitation in the same position [24] . The effect of the rotation of the wheelset has been neglected, which is a reasonable assumption for the vehicle speeds studied in this work and for impact simulations. This hypothesis is confirmed in Ref. [25] . The influence of the missing material on the circumferential variation in mass moment of inertia due to a discrete tread defect (in the order of 20 g for a 2 mm ellipsoid-shaped defect with semiaxis lengths of 30 mm and 24 mm) is neglected.
The letters "W" and "R" are used in the symbols for receptances, Green's functions and displacements to clarify whether the response has been computed for the wheelset or for the rail. A different mathematical formulation than the one used for the wheelset is adopted for the track since the relative motion between the wheelset and the track has to be accounted for. In a first step, ordinary Green's functions are computed from the direct and cross-receptances at several track positions to capture the response of the track up to a sufficient distance from the contact point. For a set of receptances R R , x 0 , x n ( f ) computed at n sampling points x n away from the loading position x 0 , their respective Green's functions G R , x 0 , x n (t) can be computed as [26] : This means that many more receptances have to be extracted for the track model compared with the wheelset model. However, some simplifications can be introduced for track models that are symmetric with respect to the track centre and/or are periodic. For track models where the track properties are periodic, it is sufficient to compute receptances starting from the load application point and moving the load within one sleeper bay in just one direction. If the track model is also symmetric, the direct and cross receptances computed by applying the load on one rail are identical with those computed for the other rail.
For a given train speed v , samples from the set of ordinary Green's functions G R , x 0 , x n (t) are combined to form the discrete [27] , which account for the motion of the contact point along the rail. The procedure is described in Ref. [21] . To construct discrete moving Green's functions for N R number of samples, the relation between time increments t and space increments

Simulation of dynamic wheelset-track interaction
Let F 1 (t) be the force acting in the contact point of wheel 1 and F 2 (t) the force acting in the contact point of wheel 2. The displacement ξ W 1 (t) of wheel 1 can be computed as, see also Ref. [29] : Similarly, for the displacement ξ W 2 (t) of wheel 2: where G W 1 , 2 = G W 2 , 1 due to reciprocity. A similar procedure is adopted for the two rails, with the difference that here the moving rail Green's functions are used. If the cross-coupling between the two rails via the sleepers and ballast (through the rail pads) is accounted for, then the displacement ξ R 1 (t) of rail 1 can be expressed as: Similarly, for the displacement ξ R 2 (t) of rail 2: As the track model studied in the present work is symmetric with respect to the centre of the track, the Green's functions The numerical implementation of Eqs. (4) -(7) is described in more detail in [21] .

Contact formulation
According to the procedure described in Ref. [21] , the normal contact problem is solved by the variational method using the active set algorithm proposed by Kalker [14] . It is assumed that the wheel and rail can be approximated by elastic halfspaces in the vicinity of the contact patch. At the interface between the two bodies, a potential contact area (which is larger than the actual contact area) is defined. This area is discretised into a mesh of N rectangular elements with side lengths x in the longitudinal direction and y in the lateral direction [21] . During the simulation, the potential contact area moves along the longitudinal direction at vehicle speed v .
The contact condition for each element in the contact area is formulated as where, for element I in the contact mesh, d m Iz is the vertical distance between wheel m and rail m ( m = 1, 2) and p m Iz is the contact pressure. The elements in the potential contact area are divided into an in-contact set (where d m Iz = 0 and p m Iz ≥ 0 ) and an active set (where d m Iz ≥ 0 and p m Iz = 0 ). In each time step, the elements that are in contact and their contact pressures are determined through an iterative procedure using the inequalities in Eq. (8) . For a smooth wheel-rail contact without surface irregularities, the kinematic constraint for the I -th element in the contact patch between wheel m and rail m reads [10] : where z m R ,I and z m W ,I are the profiles of rail m and wheel m , and ξ R m and ξ W m are their respective vertical displacements calculated according to the procedure in Section 2.2 . Note that the lateral position of the wheelset is prescribed as constant throughout the simulation. ξ W m accounts both for the elastic deformations of the wheelset and the compression of the suspension/rigid body motion of the wheelset. However, concerning the elastic deformation in the contact zone, the mesh is not fine enough to give realistic local instantaneous displacements in the contact zone. The first value in the wheel Green's functions contains the instantaneous displacements in the contact point stemming both from the (unrealistic small) local deformation in the contact zone and the compression of the suspension. In the solution procedure, this first value of the Green's functions is set to zero, since this simplifies the solution procedure considerably (current wheel displacements only depend on previous known contact forces). This is an approximation.
The local instantaneous displacements in the contact zone are handled in the contact model u m Iz instead. A sketch of a track and a wheel with their respective profiles z m R ,I and z m W ,I and vertical displacements ξ R m and ξ W m is shown in Fig. 1 . The local vertical displacement u m Iz , which is the displacement difference between the rail and the wheel, is given by where I and J are elements of the potential contact area, and A IzJz is the element of the influence coefficient matrix given by Kalker [28] describing the vertical displacement at the centre of element I due to a normal unit pressure on element J .
In each time step, the total vertical contact force F m (t) acting between wheel m and rail m is computed by summing the contributions from the contact pressure values in all elements of the contact patch.
To account for discrete irregularities, such as wheel flats [10] and rail squats [12] , their geometry is added to the contact formulation. It is also possible to account for other types of deviations from the nominal geometry, such as wheel and rail roughness or wheel out-of-roundness. The above-mentioned irregularities are discretised for each element of the potential contact area and accounted for in the extended contact formulation in Eq. (12) . For the I -th element of the potential contact area between wheel m and rail m , the total contribution in terms of deviations from the nominal geometry of wheel and rail due to roughness and discrete defects are denoted by r m W ,I and r m R ,I . Then, Eq. (9) becomes d m Some implementations of this problem and their results will be presented in Chapter 5.

Calculation of stress
Based on the calculated time histories of external forces acting on the wheelset, the time-variant stress at critical locations in the wheelset is calculated in a post-processing step. A consistent approach to accomplish this is to apply a convolution integral formulation based on pre-computed stress frequency response functions (stress over applied force [Pa/N]).
These stress frequency response functions s k,m describe the stress response in the frequency domain at wheelset position k due to the wheel-rail contact force applied on wheel m . The wheelset model applied in this paper, see Section 3 , includes the primary suspension modelled as two sets of one spring and one viscous damper coupled in parallel, one at each end of the wheelset axle. The position of the upper end of each primary suspension model is constrained. Thus, the contribution from the load in the primary suspension to the stress response at position k is accounted for in the frequency response function. Other external loads acting on the wheelset can be braking forces on the brake discs or power transmission loads imposed on the gearbox. In the present study, where the objective is to study the influence of wheel tread defects on axle bending stress, the additional stress contributions from braking and power transmission will be neglected.
Here, the stress frequency response functions are calculated by means of the subspace-based steady state dynamic analysis procedure in the commercial FE software. The complex-valued transfer functions, for each of the six stress components ( σ xx , σ yy , σ zz , σ xy , σ xz , σ yz ), are computed at several points on the axle where high stresses are anticipated.
The time history for a specific stress component σ i j in the studied point k is determined by convolving the time histories of the contact forces F 1 (t) and F 2 (t) with the stress Green's functions S k, 1 i j and S k, 2 i j for the contact points on wheels 1 and 2 The Based on the six stress time histories calculated at point k , the time history of the stress invariants in the studied nodes can be determined.

Wheelset model
The generated FE model is based on the design of a powered wheelset mounted under one end of the bi-directional X40 Swedish regional passenger train. To simplify the model, minor chamfers are not accounted for. Moreover, the geometry of the gearbox is approximated as a hollow cylinder with mass 200 kg under the assumption that 40% of the whole gearbox weight is borne by the axle. The wheels are modelled with a worn diameter of 863.5 mm (nominal diameter is 920 mm) and the web mounted brake discs are included. The wheelset model is characterised by a non-symmetric design since the gearbox is located closer to one of the wheels, see Fig. 2 .
The wheelset model was obtained by assembling six separate parts: the axle, two wheels, two pairs of brake discs and the gearbox. The axle boxes are modelled as reference points with an inertial mass of 208 kg located at the middle of the journals. The primary suspension is connected to the axle via the axle box reference points. The parts were tied at their contact surfaces, thus disregarding any connector features (bolts, interference fit). The neglected influence of the interference fit between components attached on the axle causes an underestimation of the stresses at nearby fillet radii [30][31] . As discussed in Section 2.4 , torsional stresses in the axle due to traction or braking are neglected here.
The wheelset model was meshed using twenty-node quadratic brick elements referred to as C3D20 in Ref. [23] . The average element side length in the axle is 3.5 cm, but finer elements (with side lengths of 1.5 cm) are used at the fillets in the vicinity of the wheel and gearbox seats and in the vicinity of the nominal rolling circles on the wheel treads. The highest stresses are expected close to the fillets [31] because of stress concentration effects and, as a consequence, these  finer elements were required for the extraction of representative stress frequency response functions. The model has 106,600 degrees of freedom (dofs).
As the mesh was swept with respect to the axle axis of rotational symmetry, the average element side length in the wheels, brake discs and gearbox is larger than 3.5 cm. This implies that accurate stress frequency response functions may not be captured for these components unless a finer mesh is used. However, the computation of stresses outside of the axle is not within the scope of the present study.
The primary suspension is modelled by a spring and a dashpot acting on each axle box, with stiffness 2 MN/m and viscous damping 20 kNs/m in the vertical direction. As mentioned in Section 2.4 , to apply a prescribed static wheel load, the position of the upper end of each primary suspension model is constrained. Thus, the forces in the primary suspension balance the wheel-rail contact forces by vertical loads on the axle journals. The same steel properties are used in the whole assembly (density 7810 kg/m 3 , Young's modulus 210 GPa, Poisson's ratio 0.30). The mass of the assembly is 1943 kg. Material damping is modelled using Rayleigh damping [32] . The parameters describing the stiffness proportional damping (3 ·10 −8 ) and mass proportional damping (2 ·10 −4 ) were chosen such that the average values of material damping for wheel modes involving different numbers of nodal diameters agree with the modal damping ratios for a wheel model suggested in [33] . At lower frequencies, most of the damping for the wheelset is supplied by the primary suspension.
Frequency response functions, see Fig. 3 , were extracted according to the procedure described in Section 2.1 . The corresponding Green's functions are shown in Fig. 4 . To determine the positions of the highest stresses in the wheelset axle, a separate dynamic analysis of the overall stress field in the wheelset was run in the commercial FE software. Wheel-rail contact forces computed in the in-house software were used as inputs. The forces were computed for the case of one wheel with a discrete wheel tread defect (depth 2 mm and longitudinal and lateral semi-axis lengths 30 mm and 24 mm, respectively). In Fig. 2 , the positions with the highest stresses are labelled A1, B1, C1. Here, A1 is at the fillet radius near the wheel seat that is farthest away from the gearbox), B1 is at the fillet radius near the seat on the external side of the gearbox), while C1 is at the fillet radius near the wheel seat on the gearbox side. Stress frequency response functions computed at node A1 are shown in Fig. 5 . To verify the quality of the axle mesh, two alternative models were developed where the size of the elements located near node A1 was halved once or twice relative to the mesh in Fig. 2 . Based on static analysis, the calculated difference in maximum bending stresses obtained using the mesh shown in Fig. 2 and the finest mesh was below 2%, and overall a good convergence in the results was observed. This confirms that the axle mesh used in the present study has sufficient resolution at node A1.

Track model
A 3D periodic FE model of a tangent track section has been established to calculate the moving Green's functions that are used as input in the simulation of the dynamic wheelset-track interaction. The model accounts for the cross-coupling between the two rails via rail pads, sleepers and ballast. A parameterised FE model of the ballasted tangent track has been generated in the FE software using Python scripts, using the same approach as described in Ref. [34] . The mass, stiffness and damping matrices generated by the FE software are imported to a Matlab software, where the frequency response functions for the track are calculated. Only the vertical track dynamics is considered.
The modelled track section has to be long enough to satisfy two requirements. First, a sufficient number of sleeper bays has to be present so that cross-receptances can be extracted at a sufficient distance from the excitation point. This means  that moving Green's functions can be established for a longer distance, which is required for simulation of high train speeds. Secondly, the track model needs to be long enough such that reflections of vibrations from the boundaries of the modelled track are negligible at the studied central part of the track model. Ref. [35] includes a derivation of non-reflective boundary conditions that could be easily extended to a periodic track model. If non-reflective boundary conditions are used, the model does not need to be long enough to avoid the reflections of vibrations. However, as the track model used in the present work is made of beams and linear springs and dampers, the computational requirements to extract frequency response functions and compute moving Green's functions are low. As a consequence, the length of the model used in the present study was extended to 140 sleeper bays. The constant sleeper spacing is 0.60 m and the distance between the two rails is 1500 mm. Input data to the track model are listed in Table 1 . Each sleeper is modelled as a uniform Rayleigh-Timoshenko beam discretised using 20 elements of equal length. Each rail is modelled as a Rayleigh-Timoshenko beam with a discretisation of 8 elements per sleeper bay. The rail ends on the two sides of the track model have clamped boundary conditions. Fig. 10. Plot of the variation in wheel-rail contact forces due to a discrete wheel tread defect. (a) Irregularity at lateral centre of discrete wheel tread defect on wheel 1, and (b) calculated time history of vertical wheel-rail contact forces for both wheels. Distance is measured from the point where the longitudinal centre of the defect is aligned with a position on the rail directly above a sleeper. Size C defect with h = 2 mm (see Table 3 The stiffness and damping of the ballast and subgrade are modelled by non-interacting springs and viscous dampers acting in the vertical direction between each sleeper and the ground. For each spring-damper pair, the ground is modelled as a rigid reference point located exactly below the corresponding sleeper node. All the translational and rotational degrees of freedom of these reference points are locked. The total ballast stiffness and viscous damping per sleeper is k b = 104 MN/m and c b = 86 kNs/m, respectively. A two-dimensional sketch of a cross-section of the track model is shown in Fig. 6 .  Fig. 1 , non-symmetric excitation). Longitudinal centre of defect is aligned with a position on the rail directly above a sleeper. For all rail nodes in the central 40 sleeper bays, frequency response functions (magnitude and phase angle of the vertical receptances) have been extracted up to 2500 Hz with a constant spacing of 0.1 Hz, see Fig. 7 . Two examples of moving Green's functions are plotted in Fig. 8 .

Demonstration examples -Wheel tread damage
In this Section, the presented method to simulate the vertical dynamic vehicle-track interaction is demonstrated. In particular, the influence of a discrete form of wheel tread damage on the vertical wheel-rail contact forces and axle bending stresses is studied by superimposing a discrete tread defect on the nominal geometry of the wheel tread during one wheel rotation (the wheel is assumed to have a cylindrical profile in the following simulations). The mathematical formulation for a 3D discrete defect proposed in [12] has been adopted, aiming to represent a breakout of tread material due to RCF damage.
where a and b are the longitudinal and lateral semi-axis lengths and h is the depth of the defect. An example of a 3D discrete defect described by the formulation of Eq. (22) is shown in Fig. 9 . Three different combinations of a and b will be studied in the paper, see Table 2 . In all simulations, the tread defect is placed on the nominal rolling circle of the wheel such that the lateral centre of the defect is aligned with the lateral centre of the rail. Further, the longitudinal centre of the tread defect is aligned with a position on the rail that is directly above one sleeper. All simulations feature the wheelset model in Section 3 with axle load 20 tonnes running over the tangent track model in Section 4 at constant speed 100 km/h (except for the simulations in Section 5.5 ).

Time history of wheel-rail contact force and axle bending stress
In this sub-chapter, a size C defect with depth h = 2 mm is present on one of the wheels (non-symmetric excitation). The defect is placed on wheel 1, cf. Fig. 2 . The applied Green's functions represent the dynamics of the wheelset and track up to 2500 Hz. Fig. 10 shows the calculated time history of the wheel-rail contact forces for the two wheels for a speed of 100 km/h. A preloading phase has been applied in the computation of the wheel-rail contact force to allow for the primary suspension to reach an equilibrium, see Ref. [11] . As the wheel-rail contact reaches the leading edge of the tread defect, the wheel is subjected to an unloading phase. In this phase, the wheel is moving downwards and the rail upwards to compensate for the removed wheel material. For this combination of train speed and defect size, the contact force on the wheel with the tread defect is reduced to zero over a short time interval, which means that loss of contact between the wheel and the rail has occurred (minimum contact force = 0 kN). At the far end of the defect, the contact force increases abruptly due to a reversal in the vertical wheel trajectory (not shown here). For the damaged wheel, the maximum contact force is 358 kN. After the occurrence of the maximum force, a transient characterised by minor oscillations in magnitude is observed, before the contact force returns to its original level. Note that the wheel without the defect is only subjected to a minor variation in contact force, indicating that the cross-coupling between the contact points on the two wheels is relatively weak.
To calculate axle bending stresses according to the procedure described in Section 2.5, the time histories of the two contact forces in Fig. 10 have been convoluted with their respective stress Green's functions computed for node A1. The resulting time history of the maximum principal stress is plotted in Fig. 11 . Here, the maximum and minimum stresses at position A1 are 93.9 MPa and 38.3 MPa, respectively. Before the defect enters the wheel-rail contact, the oscillation in axle stress is dominated by the sleeper passing frequency and by the effect of various wheelset modes, see Fig. 12 . The wheelset mode illustrated in Fig. 12 (a), which occurs at a frequency of 206 Hz, has a significant effect on the bending stresses in the axle, cf. Fig. 4 . In Fig. 9 , the oscillations due to this bending mode are particularly evident after the wheel-rail impact. The lower-magnitude oscillations at higher frequency are due to a wheelset mode at 1920 Hz, see Fig. 12 (b).

Convergence study
As described in Chapters 3 and 4, frequency response functions and Green's functions have been calculated accounting for frequencies up to 2500 Hz. To investigate if this frequency range is sufficient, a convergence study has been performed. Different sets of frequency response functions and Green's functions for the wheelset and track models were generated for maximum frequencies up to 2500 Hz at intervals of 250 Hz. The studied load case is the same as in Section 5.1 .
For each maximum frequency used in the analysis to generate the Green's functions, the calculated maximum and minimum wheel-rail contact forces are plotted in Fig. 13 (a). Similarly, stress frequency response functions were calculated for the six stress components at node A1 of the wheelset. These were used to obtain different sets of stress Green's functions with the same upper frequency limits as those used for the contact points displacements. The calculated maximum and minimum axle bending stresses at the investigated node are plotted in Fig. 13 (b). The results show a significant influence of frequency range on the maximum wheel-rail contact force and on the calculated axle bending stress (respectively around 20% and 15% if the frequency limit is increased from 250 Hz to 2500 Hz). It is concluded that variations in both wheel-rail contact forces and axle bending stresses are captured accurately if frequencies up to at least 10 0 0 Hz are accounted for in the analysis. Nevertheless, in the remaining parts of this paper, the frequency range up to 2500 Hz has been considered.
To obtain the required time resolution t = x/ v of the Green's functions for each maximum frequency in Fig. 13 , the length of the corresponding frequency response vector was extended by zero padding from the maximum frequency up to half of the required sampling frequency 0 . 5 f s = 0 . 5 / t (which is 13,900 Hz for train speed 100 km/h and contact mesh discretisation 1 mm).

Non-symmetric loading
Wheel-rail impact loads and axle bending stresses have been calculated for three different configurations of the size C defect with h = 2 mm. The calculated maximum stresses at nodes A1 -C1 are compared to those for the case without any tread defect. The four studied cases are: (1) no discrete defect on either wheel, (2) defect on wheel 1 (cf. Fig. 2 ), (3)  defect on wheel 2, and (4) defect on both wheels. The wheelset is thus subjected to a symmetric excitation in cases 1 and 4 and to a non-symmetric excitation in cases 2 and 3. The analyses were performed for axle load 20 tonnes and train speed 100 km/h. The results for all the studied cases are collected in Table 3 .
The results show that, if the excitation is symmetric, the difference in resulting maximum wheel-rail contact forces at the two contact points is lower than 4%. Moreover, the maximum force witnessed by a damaged wheel is very similar if one wheel is excited (damage on one tread) as when both wheels are excited (damage on both treads). In fact, the simulation results show an increase of only 2 kN for such cases. Overall, the effect from the non-symmetric wheelset design is not negligible, while the effect from the modelled cross-coupling in the wheelset and in the track is rather weak. However, for the non-symmetric damage positions, it is observed that the maximum wheel-rail contact force increases slightly for the non-damaged wheel. This is due to the cross-coupling in the wheelset and in the track.
Further, it can be observed that the wheel seat fillet (node A1) is subjected to higher stresses than the fillet at the gearbox seat (node B1), irrespectively of the position of the tread damage. The stress increase at point B1 is between 11% and 28% for all cases of tread damage compared to the no damage case. For point A1, it is between 25% and 39%, with the highest increase found when the damage is on both wheels. If no damage is present on either wheel, the fillet between the Fig. 17. Influence of depth of tread defect on maximum principal stress generated at the wheel seat fillet (node A1 in Fig. 2 ). Same load case as in Fig. 14 . gearbox and the nearby wheel (node C1) is subjected to a lower maximum stress than the maximum value for the other two fillets. However, for symmetric excitation with damage on both wheels, the maximum stress at node C1 is significantly higher than at the other studied nodes. Compared to the reference case without tread defects, node C1 is for this case subjected to the largest relative increase in stress, with an increase by 81%. The effect from the modelled cross-coupling is significant for the axle stresses.
It is concluded that the studied combinations of tread damage lead to significant wheel-rail impact loads (high dynamic load magnification factors relative to the static wheel load). The corresponding relative increase in maximum axle bending stresses is not as pronounced. It is argued that this is due to the damping and inertia of the wheelset assembly, which lead to an attenuation of the high-frequency contributions to the axle stresses.

Size of wheel tread defect
Another parameter study has been performed to study the influence of tread defect size on the wheel-rail impact load and axle bending stress. The three combinations of a and b listed in Table 2 are studied. For each given shape of the defect, the depth was increased from 0 mm to 3 mm at regular intervals of 0.05 mm.
For axle load 20 tonnes and train speed 100 km/h, the maximum and minimum vertical wheel-rail contact forces due to a defect placed on the nominal rolling circle of wheel 1 are plotted in Figs. 14 (a) and 14(b). Fig. 14 (a) shows a substantial increase in impact load with increasing depth of the defect h , in particular for sizes B and C. For some values of h , it is observed that a smaller defect (in terms of a and b ) may cause a higher impact load than a larger defect, cf. the size A and C defects for h < 0.13 mm or the size B and C defects for h < 0.44 mm. This is due to a more pronounced change in the vertical trajectory followed by the centre of the wheel that may occur for smaller defects, see Fig. 15 .
In Fig. 14 , it is observed that the impact load converges towards a constant value as the depth of the defect exceeds a certain value, which depends on the combination of a and b and the nominal wheel-rail contact geometry. This is because the vertical trajectory followed by the centre of the wheel will no longer be influenced by an increasing defect depth when the wheel-rail contact is instead supported by the contact around the edges of the defect (above a certain depth there will not be contact at the bottom of the defect), see Fig. 16 . Fig. 16 also shows that different sizes of the defect may lead to very different distributions of the contact pressure even if the depth of the damage is the same. The size B and C defects lead to loss of wheel-rail contact for h larger than about 0.65 mm.
In a post-processing step, the time histories of the calculated contact forces have been used to compute time histories of axle bending stresses using the convolution procedure described in Section 2.4 . This procedure has been applied for the cases shown in Fig. 14 using the stress Green's functions for node A1. The maximum principal stresses for these cases, when evaluated over a travelling distance of three sleeper bays, are plotted in Fig. 17 . Since only the vertical component of the wheel-rail contact force is accounted for, the value of the maximum principal stress is dominated by the value of the axle bending stress.
For the size A defect in Fig. 17 , there is no visible effect on the axle bending stress for any of the studied depths. In this case, the time history of the axle stress is dominated by the parametric excitation from the sleeper passing. For the size B and C defects with increasing depth, the increase in stress becomes significant as the influence of the wheel-rail impact load starts to exceed that of the parametric excitation. For the size B defect, the increase in maximum stresses is in the order of 12 MPa. For the size C defect, the increase in maximum principal stress is more significant, reaching about 25 MPa.

Train speed
For non-symmetric loading with the three defects listed in Table 2 , the influence of train speed on the wheel-rail impact load has been studied. Two depths 0.3 mm and 2.0 mm are compared, and the defect is positioned on wheel 1, cf. Fig. 2 . Train speed is varied between 40 km/h and 220 km/h at intervals of 10 km/h. The calculated maximum and minimum contact forces are shown in Fig. 18 .
As expected, the maximum contact force is generally increasing with increasing train speed and with increasing depth of the defect. For depth 0.3 mm, it is observed in Fig. 18 (a) that the size B defect generates higher impact loads than the larger size C defect for train speeds lower than 150 km/h. On the other hand, for depth 2.0 mm, the size C defect generates higher impact loads than the size B defect for all speeds.
In a similar way as in Section 5.4 , the axle bending stress at node A1 has been calculated for the cases in Fig. 18 . The results are shown in Fig. 19 . It is observed that maximum stress does not have a similar speed dependence as the maximum contact force.
For h = 0.3 mm, the influence of speed is more pronounced above 140 km/h, in particular for the size B and C defects. For h = 2.0 mm, the size A defect shows basically no change as compared to the 0.3 mm case. A more pronounced variation with speed is found for defects of sizes B and C. For an increase from 40 km/h to 70 km/h, the maximum principal stress is increased by about 12 MPa. Local maxima in the calculated axle stresses are observed around speeds of 90 km/h, 120 km/h and 180 km/h.

Conclusions
A method to simulate the vertical dynamic interaction between a wheelset and track, capable of accounting for generic distributions and shapes of wheel tread damage or out-of-round wheels, has been presented. The dynamic behaviour of the wheelset (two wheels, axle and primary suspension, and possibly including brake discs, gearbox assembly and axle boxes) and track (two rails discretely supported by rail pads and sleepers on ballast) is assessed using 3D finite element models. The coupling between the two wheel-rail contacts (one on each wheel) via the wheelset axle and via the sleepers is considered. This means that two contact problems are simulated simultaneously, one for each nominal wheel-rail contact point in the same wheelset. The dynamics of the wheelset and track are treated using a Green's functions approach, while the wheel-rail normal contact is solved using Kalker's variational method. Wheelset designs that are non-symmetric with respect to the centre of the axle, non-symmetric distributions of tread damage on the two wheels (or irregularities on the two rails), as well as track support conditions that are non-symmetric with respect to the centre of the track (not studied in this paper) can be investigated.
Further, based on a convolution integral formulation, a procedure for computing stresses in the wheelset in a postprocessing step has been introduced. For various combinations of train speed and size of tread damage, the method has been demonstrated by calculating wheel-rail impact loads and maximum axle bending stresses at different locations in the axle. The parametric study showed that the increase in wheel-rail impact loads depends on the combination of the defect shape (semi-axis lengths) and depth. An increase in train speed will lead to a substantial increase in wheel-rail impact load, while the influence of tread damage distribution and train speed on axle bending stresses is not as pronounced. For the type of discrete tread defect implemented in this paper, only defects of larger size (sizes B and C in Table 2 ) tend to have a significant influence on axle stresses, as the dynamic influence of a small defect is comparable to the influence of the parametric excitation due to the sleeper passing. In addition, the axle bending stresses will be influenced by irregularities in track geometry (not studied here).
The stress calculation procedure can be extended in a straightforward way by also accounting for stresses introduced by the interference fit between wheel and axle or by torsional loads from the gearbox transmission. In this way, stress time histories at critical locations can be determined and fatigue stress spectra can be established for various axle designs and vehicle running conditions on tangent track. Since the simulation model is based on pre-calculated Green's functions, it allows for versatile and computationally costefficient numerical predictions. This is an attractive feature for the understanding of the consequences of wheel tread damage, and for the optimisation of maintenance planning.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.