Comparisons between beam and continuum models for modelling wheel-rail impact at a singular rail surface defect

A singular rail or wheel surface irregularity, such as a squat, insulation joint or wheel ﬂat, can cause large wheel- rail impact force. Both the magnitude and frequency content of the impact force need to be correctly modelled because they are closely related to the formation, deterioration and detection of such irregularities. In this paper, we compare two types of commonly used wheel-track interaction models for wheel-rail impact problems, i.e., a beam and a continuum ﬁnite element model. We ﬁrst reveal the diﬀerences between the impact forces predicted by the two models due to a typical rail squat using a time-frequency analysis. Subsequently, we identify the causes for the diﬀerences by evaluating the eﬀects of diﬀerent model assumptions, as well as diﬀerent model parameters. Results show that the impact force consists of a forced vibration peak M1 followed by free vibration related oscillations with three dominant frequencies: f 1 (340 Hz), f 2 (890 Hz) and f 3 (1120 Hz). Compared with the continuum model, the beam model with a Hertzian contact spring overestimates the M1 peak force. The discrepancy can be reduced by using a Winkler bedding contact model. For the track model, the beam model is comparable to the continuum model up to about 800 Hz, beyond which the track damping starts to deviate. As a result, above 500 Hz, the contact forces dominate at f 2 for the beam while at f 3 for the continuum model. Finally, we show that the continuum model is more accurate than the beam model by comparing to ﬁeld observations. The eﬀects of stress wave propagation on the diﬀerences are also discussed.


Introduction
Wheel-rail vertical impact usually occurs at short wavelength defects (e.g., squats, poor welds, wheel flats, short-pitch corrugations) or structure discontinuities (e.g., insulated joints, crossings). The appropriate modelling of structural flexibility of the wheel-track system is of major concern for wheel-rail impact problems due to their high frequency nature [ 1 , 2 ]. Based on the assumptions adopted, two types of wheel-track interaction models are commonly used for wheel-rail impact problems. The first type is referred to as the beam model in this paper, in which rails are modelled using the Euler-Bernoulli or Timoshenko beams [3][4][5] . The rail can be continuously supported on the elastic or Winkler foundations [6][7][8] , a layer of sleeper beam [ 9 , 10 ], discretely supported on sleepers [11][12][13] or on substructures, such as bridges [14] or subgrades [15] . The second type drops the assumptions made in the beam theories and treats the components of the wheel-track system (i.e. the sleeper, rail, wheel or even railpads) as continua. This type of models is referred to as the continuum model in this paper. Continuum mod-length of squats is about 20-40 mm, e.g., in the Dutch railway network [25] . With such a geometry, the contact force is characterized by the continuous short wavelength peaks in the time domain [26] . In the frequency domain, two major frequencies at around 300 Hz and 1000 Hz have been identified in the simulations and verified by the axle-box acceleration (ABA) measurements [27] . Both the magnitude and frequency content of the impact force need to be correctly modelled because they help better understand the formation [ 25 , 28 ], development [29] and detection [30] of squats.
Different wheel-track interaction models have been compared through benchmark tests [31][32][33] . For example, six different beam models were benchmarked against each other in [32] . The P1 forces at a wheelflat were found to be ranging from 2 to 7 times the static load. The track vibrations reproduced by the six models were not consistent between each other or compared to the field measurement [33] . Although the discrepancies can be pinpointed by the benchmark tests, it is nonetheless difficult to deduce exactly which assumption contributes to an identified discrepancy, because there are usually more than one different assumptions between two benchmark participants.
Alternatively, comparative studies have been carried out by changing one assumption at a time per each track component to identify their effects on the dynamic characteristics. For rail models, the Euler-Bernoulli beam is considered adequate below about 500Hz due to the neglecting of the shear deformation [ 1 , 34 ], whereas the Timoshenko beam is accurate up to about 1500 Hz [35] , beyond which the assumption of rigid cross section in the beam theories becomes invalid and a continuum model should be adopted. Moreover, the Euler beam theory tends to overestimate the P1 force compared to the Timoshenko beam theory as well as the measurement [36] . Different railpad/fastening models have been compared in terms of the effect of the support length and size for beam models [ 37 , 38 ] and for continuum models [ 39 , 40 ]. The introduction of sleepers as discrete supports has two major effects. On one hand, the pin-pin resonance can be correctly modelled using the discretely supported track [9] . On the other hand, the sleeper layer could effectively reduce the P1 force compared to the beam on elastic foundation model [41] . In the case of soft railpads, sleepers can be modelled as rigid masses [ 23 , 42 ]. As railpad stiffness increases, the coupling between the rail and sleeper becomes stronger. As a result, flexible sleeper models should be adopted to model extra track resonances due to sleeper bending [ 23 , 42 , 43 ].
Except for the track models, different wheelset and contact models have also been compared. For beam models, the wheelset can be modelled as a rigid mass or a flexible body through modal superposition [ 44 , 45 ]. Flexible wheelset models can slightly reduce the impact force at wheelflats compared to rigid models [44] . In addition, differences can be observed in the frequency domain of the contact force, limited to narrow frequency bands corresponding to the eigenfrequencies of the flexible wheelset [45] . Different contact models are compared in [46][47][48][49] . However, all these contact models are quasi-static, which do not consider local inertia effects within or in the vicinity of the contact patch. In contrast, the continuum model takes in to account local structural vibrations in the vicinity of the contact patch and stress wave propagations in the wheel and rail [26] . For example, in [50] , wheel-rail contact induced Rayleigh waves have been reproduced by a continuum model.
In the literatures, despite the extensive comparisons between different beam models, comparisons between beam and continuum models can only be found in a few cases, such as in the unloaded condition [ 23 , 42 ] or under parametric excitations [51] , but not for wheel-rail impact problems. Even for beam models, the comparisons are almost all based on wheel-rail impacts at wheelflats or joints, which are generally larger geometrical irregularities than squats. As a result, the focus was to compare the magnitudes of the contact force, e.g., the P1/P2 forces, whereas the frequency or time-frequency characteristics of the impact force have not been fully explored. Furthermore, the effects of the local inertia in the contact patch and the stress wave propagation in solids on the wheel-rail impact force have not been studied. This paper aims to compare a beam and a continuum model for simulating the impact forces induced by squats. To this end, two FE models, meshed with beam and solid elements, respectively, have been developed (Section 2). We first identify the major characteristics and discrepancies of the wheel-rail impact force simulated by the two models at a typical squat, in the time, frequency and time-frequency domain (Section 3.1). Subsequently, we examine the effects of the model assumptions made in the wheel, contact and track models separately (Sections 3.2), as well as the effects of different model parameters (Section 3.3). In Section 4, a frequency domain model is adopted to investigate the coupled effects of the three components. Finally, we compare the simulated results with field observations and discuss in detail the influence of wave propagation (Section 5).

Beam model
The track is represented by a two-layer discretely supported model, see Fig. 1 (a). The rails and sleepers are meshed with the Timoshenko beam element. At each node, only the vertical and in-plane rotational degrees of freedom are considered. By a convergence analysis, the optimal mesh sizes for the rail and sleeper are determined, with 24 elements per sleeper span for the rail and 20 elements per sleeper. Ballast and railpads are modelled as discrete spring-damper pairs.
Two track model options are made: the full track model and the half track model. The latter is the former halved along the track central line, where the symmetrical boundary conditions are applied. The wheelset is simplified as a rigid body in the full track model and as a rigid mass in the half track model. The bogie and car body are simplified as static loads applied vertically on the wheel.
Only the vertical wheel-rail contact is considered with two options. The first one is a non-linear Hertzian spring model, for which the halfspace assumption applies. The contact force is calculated as are the vertical coordinates of wheel, rail and defect geometry, respectively, and C H is the Hertzian coefficient and can be approximated as where E and are the Young's modulus and Poisson's ratio of the wheel and rail (assumed equal for both materials), and R is the radius of the rail head in the lateral direction. The defect geometry is analytically defined as a cosine function: where D, L and x 0 are the depth, the length and the starting location of the defect. The second contact model is a Winkler bedding model [60] which uses multiple springs in the longitudinal direction, see Fig. 1 The wheel-track interaction model is solved in the time domain using the Newmark integration with a fixed time step length of 4e-5 s. To ensure the convergence for the contact force, the Newton-Raphson iteration is adopted within each time step. The solution process is implemented in Matlab.

Continuum model
The second model is a detailed three-dimensional (3D) FE model meshed with continuum elements, as shown in Fig. 2 . To reduce the computational cost, only halves of the track and wheelset with the symmetrical boundary condition are considered. The bogie and car body are simplified as sprung masses. The railpads and ballast are modelled as discrete spring-damper pairs. For each railpad, a grid of 3 by 4 springdamper pairs is adopted between the rail and sleeper, see Fig. 2 . The rolling of the wheelset is enabled by applying a torque on the axle. This model is developed and solved with the commercial FE program ANSYS/LS-Dyna. An 8-node hexahedra element with reduced integration points is adopted due to its computational efficiency. The length of both track models is 12 m (i.e. 20 spans), which has been shown to be long enough to get rid of the boundary effect for this problem [26] . The wheel-rail contact is solved using a surface-to-surface contact algorithm incorporated in LS-Dyna. This algorithm with the mesh size of 1 mm in the contact patch has been shown to yield satisfactory results in comparison to the Hertz or CONTACT [52] model in the quasi-static wheel-rail contact case [53] . The longitudinal geometry of the irregularity is defined as in Equation (3) , while its lateral geometry remains the same as the lateral profile of the rail head, see Fig. 2 .
The dynamic wheel-track interaction is solved using an explicit time integration scheme (a central difference method) with a time step length of 3e-8 s, which is small enough to ensure both the stability and convergence. For the details of the modelling and solution procedure, the readers are referred to [54] .

Characterization of wheel-rail impact force
We first consider the wheel-rail impact at a typical squat in the Dutch railway network as a reference case. The model parameters are listed in Table 1 . These parameters are taken from [ 27 , 55 ], which have been validated by field hammer tests [55] and ABA measurement [27] . We model the fastening system as spring-damper pairs: railpads sustain compression and clamps and bolts sustain tension. This is a widely accepted simplification in railway track models. As in the loaded condition the fastening system is in compression, we take the railpad stiffness value as the spring stiffness. The sleeper is a prestressed concrete sleeper of type NS90. We assume a uniform cross section for the sleeper with equivalent cross sectional properties. The squat geometry is defined using Equation (3) with D = 0.2 mm, L = 30 mm and x 0 = 6.5 m. This geometry is chosen because a typical squat in the Netherlands is between 20 mm ~40 mm long [25] and less than 0.4 mm deep [56] . Besides, the squat is located near a sleeper support at 6.6 m, which is also typical in the Dutch railway network [25] . For the convenience of comparing to measurements (Section 5), the velocity of the wheel is set to 30 m/s, which is the same velocity at which the ABA was measured in [27] .    Fig. 3 (a)) consist of two stages. From 6.5 m to 6.53 m (indicated by the vertical dashed line), the wheel-track system is in a forced vibration stage, excited by the defect geometry. The contact forces show a local dip (D1) followed by a local maxima (M1). The M1 magnitudes are 2.1 and 1.4 times the static load for the beam and continuum model, respectively. From 6.53 m on, the system is in a free vibration stage. The contact forces oscillate in two major wavelengths, i.e. a shorter wave D2M2D3 and a longer wave D3D4. In the frequency domain, there are three major characteristic frequencies (see f 1 , f 2 and f 3 in Fig. 3 (b)). The magnitudes calculated by the continuum model are lower at f 1 and f 2 than the beam model. In addition, the dominant frequency above 500 Hz is at f 2 = 890 Hz and f 3 = 1130 Hz for the beam and continuum model, respectively.
The wheel-track system is not a constant but a time-variant system as the wheel position changes. Consequently, the characteristic frequencies should also be changing with the wheel position. To illustrate this, we apply the synchrosqueezed wavelet transform [57] to the contact force and plot the synchrosqueezed wavelet power spectrum (SWPS) in Fig. 3 (c)~(f). Three frequency bands are clearly visible in the plots. There are a number of differences between the SWPS calculated by the two models. The major difference lies between about 6.53 m and 6.6 m, where the largest power concentrates at f 2 and f 3 for the beam and continuum model, respectively. The two dominant frequencies correspond to the two wavelengths of D2M2D3 in Fig. 3 (a), i.e., approximately 38 mm and 27 mm for the beam and continuum model, respectively. Besides, the frequency change of f 3 is less abrupt for the continuum model than the beam model, as indicated with the white boxes in Fig. 3 (e) and (f). This is mainly because of the modelling of railpad as a grid of 3 by 4 spring-damper pairs in the continuum model as opposed to a single spring-damper pair in the beam model (see Fig. 2 and Fig. 3 (g)). This unphysical abrupt change of stiffness due to the singlepoint supported railpad model has also been observed in Timoshenko beam models [ 58 , 59 ] under parametric excitations. After 6.6 m, the two models show similar results with only f 1 and f 3 presented. The difference is that the contact force decays more quickly for the continuum model than the beam model, resulting a lower power in the SWPS after 6.6 m.     The simulation results by the two models shown in Fig. 3 may be influenced by two factors: one is the model assumptions and the other is the model parameters . To make a more comprehensive comparison of the two models, we investigate the effects of model assumptions and parameters in Section 3.2 and 3.3, respectively.

Effects of model assumptions
Compared to the continuum model, more assumptions are made in the beam model in terms of the wheel, contact and track model. In this section, some key assumptions are varied, see Table 2 , to investigate their effects on the simulation results. While varying model assumptions, we use the same set of model parameters (see Table 1 ) for the analysis in this section.
For the wheelset model, a rigid mass is assumed for the beam model used in Section 3.1 (Model A1 in Table 2 ). Compared to a flexible wheelset model, this assumption might result in an overestimation of the impact force for wheel flats [44] , as well as different frequency contents of the contact force at corrugations [45] . Therefore, in Section 3.2.1, the effects of wheelset flexibility on the impact force at squats are examined.
For the contact model, the half space assumption is made in the Hertzian spring model. However, the typical length of a wheel-rail contact patch in the longitudinal direction is approximately 15 mm, which is comparable to the typical length of squats (20 ~40 mm). This means the half space assumption may no longer be valid in the longitudinal direction. To account for local geometric variations within the contact area in the longitudinal direction, a two-dimensional contact model with multiple independent springs is considered (Model A3 in Table 2 ), referred to as the Winkler bedding model [60] . This model is compared to the Hertzian spring model and the 3D contact model in Section 3.2.2. For the track model, the main differences between the beam and continuum model are in the assumptions made for the rail and sleeper. It is conventionally believed that the Timoshenko beam is accurate up to about 1500 Hz, due to the assumption of a rigid cross section. Inherent from the different rail and sleeper models are the different fastening models. The rail is assumed to be supported at a single point in the beam model, whereas in reality the support is over an area, which is more realistically modelled in the continuum model, as shown in Fig. 2 . The effects of these assumptions made in the track models are investigated in Section 3.2.3.

Wheelset flexibility
The natural frequencies and mode shapes up to 2000 Hz of the 3D FE wheelset ( Fig. 2 ) are obtained by eigen-analysis. The modal superposition method is then used to represent the flexible wheelset, of which the total DOF of the wheelset is reduced to 27 modal coordinates. In such a way, the flexible wheelset model is incorporated in the beam model and the wheel-track interaction is solved in the time domain in the same way as the rigid wheelset model.
The contact forces calculated by the flexible and rigid wheelset model are compared in Fig. 4. The contact force magnitudes of the first two peaks obtained by the flexible wheelset model are slightly smaller than those by the rigid wheelset model ( Fig. 4 (a)). The PSD of the contact force calculated by the flexible wheelset model has three troughs at three wheel resonances compared to the rigid wheel model ( Fig. 4 (b), indicated by the dotted-line boxes). The findings agree with previous studies, see for example [44] for the contact force magnitudes and [45] for the frequency content. In general, the influence of wheelset flexibility on the contact force can be neglected. We therefore assume a rigid wheelset for the rest of the paper.

Contact models
In this section, we compare two contact models adopted in the beam model, i.e. the Hertzian contact spring (Model 1A in Table 2 ) and Winkler bedding model (Model 1C), as well as the 3D FE contact model (Model B). Results are shown in Fig. 5 . By changing from the Hertzian spring to the Winkler bedding, the M1 magnitudes, as well as the FFT magnitudes at f 1 , f 2 and f 3 , are reduced. Such decreases are mainly because that the contact filter effect [60] in the longitudinal direction can be taken into account with the Winkler bedding model as opposed to the Hertzian spring model, while it is automatically considered in the 3D FE contact model. Fig. 6 compares the contact solution of the three models in more detail. When the wheel center is on the descending or ascending edge of the defect, the contact patch centers of the 3D FE model and the Winkler bedding model do not coincide with the wheel center. In particular, when the wheel center is at the lowest point of the defect, i.e. at 6.515 m, the wheel is in contact with the rail on both the descending and ascending edge of the defect, resulting in two contact patches. In contrast, it is always single point contact and the contact point is right underneath the wheel center for the Hertzian spring model. The contact filter effect does not change the frequency content; the f 1 , f 2 and f 3 are the same for the Hertzian spring and the Winkler bedding. Additionally, with the Winkler bedding model, the dominant frequency above 500 Hz is still at f 2 , instead of at f 3 as with the 3D FE model.

Track models
To exclude the influence of the wheelset and contact models, we calculate the point receptances of the track at two locations, at 6.5 m and 6.6 m along the longitudinal direction, as shown in Fig. 7 . Fig. 8 compares the receptances of the two track models at the two locations. Three major track resonances (TRs) are observed for both models in Fig. 8 (a). TR1 and TR2 are the full track resonance and the rail resonance, respectively [55] . The peaks at around 1000 Hz are the pin-pin resonances. The receptances are comparable up to the anti-resonance (at about 800 Hz) before the pin-pin resonance. Compared with the beam model, the receptance magnitude of the continuum model at this anti-resonance is larger (i.e. less deeper). In addition, the phase change predicted by the continuum model is smaller at this anti-resonance than that by the beam model; see the positive phase change at around 800 Hz in Fig. 8 (c). Another difference is that the receptance of the beam model at x = 6.5 m shows a more distinct peak at the pin-pin resonance. These observations indicate that the damping of the continuum model is larger than that of the beam model before and around the pin-pin resonance. After the pinpin resonance, however, the beam model predicts smaller changes of both the receptance magnitude and phase at the anti-resonance around 1300 Hz. This means after the pin-pin resonance, there is larger damping in the beam model than in the continuum model. In the track receptances at 6.6 m ( Fig. 8 (b)), the pin-pin resonance completely vanishes for the beam model, whereas for the continuum model, there is still a small peak at around 1000 Hz.
To compare the damping properties of the two models, the modal properties of the two models are identified using the least square rational fraction (LSRF) method [61] based on the receptances shown in Fig. 8 (a). Fig. 9 shows the change of damping ratios with frequencies.
Around the pin-pin resonance (between about 800 Hz and 1100 Hz), the damping ratio of the continuum model is higher than the beam model. This is caused by the different fastening models. The rail in the beam model is supported with a single spring-damper pair, whereas in the continuum model the rail is supported by a grid of spring-damper pairs over an area; see Fig. 2 and Fig. 7 . After the pin-pin resonance, there is an abrupt increase of damping for the beam model, resulting in a larger damping ratio between about 1100 Hz and 1800 Hz. This is likely to be caused by the stress wave propagations in the rail, as will be discussed

Effects of model parameters
For the analysis in Sections 3.2, the same set of model parameters were used. In this section, two model parameters, i.e., the wheel speed and defect geometry, are varied to evaluate their effects on the comparison of different models. These two parameters are chosen because they have been shown to have significant influence on the wheel-rail impact at wheel flats [47] and squats [ 26 , 56 ]. Besides, as has been shown in Section 3.2.2, the impact force is also sensitive to contact models. We therefore combine the change of model parameters with different contact models, i.e., Model A1, A2 and B in Table 2 , in the subsequent analysis. Other model parameters, such as the vehicle and track parameters, are kept the same as in the reference case. Fig. 10 shows the impact forces calculated with different wheel speeds ranging from 10 m/s to 40 m/s. In the space domain, the first peak force (M1) calculated by the Winkler bedding model are almost identical to those by the continuum model for all the speeds.  The Hertzian spring model, however, always yields larger M1 force. In the frequency domain, as the wheel speed increases, high frequency components become more evident for all the three models. At high speeds, e.g., 30 m/s and 40 m/s, the FFT magnitudes become dominant at f 2 for the beam model and at f 3 for the continuum model. In general, the discrepancies between the Winkler bedding model and the continuum model become smaller as the wheel speed decreases. Fig. 11 compares the impact forces calculated by different models for different squat geometries. It can be seen that the impact force is more sensitive to the length than the depth of squats. In the space domain, the two peak forces M1 and M2 merge into one peak as the length increases. In the frequency domain, the FFT magnitudes at high frequencies (e.g., at f 2 and f 3 ) gradually diminish as the defect length increases. The FFT magnitudes calculated by the beam models are always  larger at f 2 than those at f 3 , whereas for the continuum model, the f 2 magnitude only becomes dominant for longer defects, e.g., longer than 50 mm.

Defect geometry
In both the space and frequency domain, the discrepancies between different models become smaller as the defect length increases or as the defect depth decreases. In general, beam models, especially with the Hertzian spring, are only comparable to the continuum model for longer defects, e.g., in the current case at least longer than 50 mm. It should be noted that the maximum defect depth considered here is 0.3 mm, which is typical for squats while might be small for wheel flats. Further investigations are needed as to whether the conclusions made here are still valid for larger defects.
In both Fig. 10 and Fig. 11 , despite the change of the FFT magnitudes, the characteristic frequencies of f 1 , f 2 and f 3 remain the same for each model. This suggests that f 1 , f 2 and f 3 represent certain resonances of the coupled wheel-track system. The origin of these resonances will be further investigated in Section 4.

Coupled dynamics of the wheel-track system
We have analysed the influences of the wheelset, contact and track models separately in Section 3. The three characteristic frequencies of the impact force, i.e., f 1 , f 2 and f 3, are not sensitive to the change of wheelset and contact models. For different track models, we only compared their effects on the track receptance (see Fig. 8 ); how the track receptance is correlated with the impact force still remains unclear. In this section, we use a frequency domain model to combine the effects of the three components. In such a way, the origin of the characteristic frequencies, as well as the contribution of each model component to the characteristic frequencies are clarified.

Frequency domain model
The defect geometry Z irr ( x ) defined in Equation (3) can be transferred into the time domain by dividing Z irr ( x ) by the velocity of the wheel that passes over the defect. Then the time domain defect geometry can be further transferred into the frequency domain using the Fast Fourier Transform (FFT). The Fourier coefficient at different frequencies f is denoted as X ( ), where = 2 f is the circular frequency. Fig. 12 shows the geometries of two short-wave rail surface defects with different lengths. In the frequency domain, the Fourier coefficients are flat up until a certain cut-off frequency and the larger defect has a lower cut-off frequency.
The contact force Y F ( ) can then be solved in the frequency domain assuming a linear time-invariant system as where H F ( ) is the transfer function that represents the characteristics of the wheel-track system and can be formulated as [ 9 , 62 ], where w , t , c are the point receptance of the wheel, track and contact spring at the contact point, respectively. It should be noted that Equation (5) is derived under the condition that only half of the track is considered.
To employ this frequency domain model for the analysis, we need to make some assumptions. First, a linearized contact spring stiffness K H is assumed. To account for the nonlinear behaviour of the Hertzian contact spring, i.e., the change of contact stiffness with contact force,we adopt a range of K H between 4 × 10 8 N/m and 1.2 × 10 9 N/m, according to the range of the contact force, see, e.g., Fig. 3 (a). The receptance of the contact stiffness can be calculated as c ( ) = 1/ K H [ 9 , 62 ]. Second, as the effect of wheelset flexibility is negligible (see Section 3.2.1), the wheel is considered as a rigid mass of M w = 900 kg (approximately half of a motorized wheelset). Thus the wheel receptance can be calculated as w ( ) = 1/ 2 M w [ 9 , 62 ]. Third, to represent the non-stationary feature of the wheel-track system as shown in Fig. 3 (c)~(f), we adopt the concept of ' frozen configuration ' [63] of the wheel-track system for each wheel position x . More specifically, we assume that at each wheel position x , the system is ' frozen ' (made stationary) with the track receptance denoted as t ( , x ). Examples of t ( , x ) at x = 6.5 m and 6.6 m are shown in Fig. 8 . Thus, the ' frozen ' transfer function H F ( , x ) and the force spectrum Y F ( , x ) can be calculated according to Equation (5) and (4) , respectively. In such a way, the non-stationary wheel-track system changing with the wheel position can be represented by a sequence of ' frozen ' stationary systems.

Correlating track receptance to impact force
We first calculate the frozen configuration response at x = 6.5 m. Fig. 13 shows the receptances ( ) (top row) and force spectrums Y F ( ) (bottom row) calculated at x = 6.5 m for the continuum model (first column) and the beam model (second column). From the frequency domain model, we can see that f 1 , f 2 and f 3 are due to the coupling between the different components of the wheel-track system. At these frequencies, the wheel receptance is much smaller than the track and contact receptance. Therefore, the characteristic frequencies of the contact force can be approximated by the intersection points of the track receptance curves with the contact stiffness line, neglecting the effect of wheel mass. More specifically, the contact stiffness line intersects with the receptance curve to the right side (the mass dominated part) of the three track resonances, i.e., the TR1, TR2 and pin-pin resonance, resulting in f 1 , f 2 and f 3 , respectively. Hence, f i ( i = 1, 2, 3) can be seen as the resonance frequency of a single-degree-of-freedom system, as shown in Fig. 14 , with the equivalent mass and the contact stiffness K H For comparison, the Fourier spectrums from the corresponding time domain models are also presented in Fig. 13 . The f 1 , f 2 and f 3 obtained by the time domain models agree relatively well with their counterparts obtained by the frequency domain models with K H = 1.2 × 10 9 N/m. In terms of the FFT magnitudes, the time-domain continuum model predicts lower magnitudes at f 2 and f 3 compared to its frequency domain model with K H = 1.2 × 10 9 N/m ( Fig. 13 (c)). This is because in the frequency domain model, the contact filter effect is not considered. For the beam model, in which the contact filter effect is also not considered, the FFT magnitudes agrees well between the time and frequency domain models. It is also noticed that above 500 Hz the frequency-domain model with K H = 1.2 × 10 9 N/m predicts a higher magnitude at f 3 for the continuum model whereas at f 2 for the beam model ( Fig. 13 (c) and (d)). This is in line with the time-domain predictions. This means the difference in the dominant frequencies between the two models is due to the difference in the track models. More specifically, it is due to the different damping properties of the two track models. As shown in Fig. 9 , the damping ratio above the pin-pin resonance is much larger in the beam model, which leads to the attenuation of the f 3 magnitude in Fig. 13 (d). Likewise, as the damping ratio is larger in the continuum model before the pin-pin resonance, the f 2 magnitude in Fig. 13 (c) is smaller than the f 3 magnitude. Now we examine the time-variant feature of the wheel-track system using the frequency domain model. We compare the predictions of f 1 , f 2 and f 3 by the frequency domain model with those by the time domain models, as shown in Fig. 15 . Comparing the frequency domain solutions (see the red circles in Fig. 15 ), the beam model predicts a larger variation of f 3 with the change of wheel positions than the continuum model. This is due the different support lengths in the longitudinal direction of the fastening models in the two models (see Fig. 3 (g)).
The time domain response is equal to its ' frozen ' part plus a term representing the dynamic effects [63] . Therefore, by using the frequency domain model as a baseline, we can compare the dynamic effects of the two time domain models. In Fig. 15 , the SWPSs are obtained by the time domain models. Unlike the frequency domain solutions, the time domain solutions of f 2 and f 3 are asymmetrical about the sleeper support at 6.6 m. For both the continuum and beam model, the largest deviations of the time domain solutions from the frequency domain solutions occur between 6.55 m and 6.6 m, predominately at f 3 . The deviations are larger for the beam model than the continuum model. This means the dynamic effects, such as due to wave propagations (to be discussed in Section 5.2), in the beam model are more pronounced than those in the continuum model, when the wheel is approaching the support.

Accuracy of the modelling results compared to field observations
Continuum models with similar approaches have been validated using both the axle-box acceleration measurement [27] and the wear patterns observed in the field [64] . We discuss how the simulation results in this paper fit with these measurements and observations.

Axle box acceleration
In [27] , the ABA signals were obtained with accelerometers mounted on the four axle boxes of a bogie. The sampling frequency was 25000 Hz, which is also the sampling rate we used in this paper to sample the simulation outputs. The measured ABA signals were low pass filtered with a cut-off frequency of 2000 Hz. In this paper, the simulation results were not filtered as they are less noisy than the measurement. However, as the frequency of interest is up to about 1200 Hz in this paper, the filtered ABA measurements are still valid for comparison with simulations. In addition, a uniform lateral profile of the squat is considered in this paper, see Fig. 2 . In comparison, the ABA signals were obtained at squats with non-uniform lateral profiles. Nevertheless, it was shown in [65] that the wavelengths or characteristic frequencies of the impact force are not influenced by the lateral profile of squats.
The ABA measured at various squats show two distinct frequency bands, with the lower frequencies between 300 ~500 Hz and the higher frequencies between 1000 ~1200 Hz [27] . In this paper, the two models yield nearly identical f 1 frequency at around 340 Hz. However, the higher frequency calculated by the beam model is at f 2 = 890 Hz and that by the continuum model is at f 3 = 1120 Hz. This means the continuum model fits better with the ABA measurements in terms of the higher frequency (between 1000 ~1200 Hz).

Wear pattern following squat
Corrugation-like wear patterns can be observed after squats in the train traffic direction [ 25 , 64 ]. Examples of such patterns are shown in Fig. 16 . The direct dynamic effect of a squat is the first wave pattern after it, which is usually shorter than 30 mm. This wave is caused by the first peak of the contact force after the defect, i.e. the D2M2D3 in Fig. 3 . In this paper, the continuum model gives a more accurate prediction of this wavelength (26 mm), while the beam model overestimates the wavelength (38 mm).

Stress wave propagation in track
One major dynamic effect that is more realistically modelled in the continuum model than in the beam model is the stress waves propagat-ing in solids. In general, there are three types of waves due to dynamic loadings, i.e. Rayleigh waves, shear waves (S-waves) and dilatational waves (P-waves). Rayleigh waves propagate near surface while S-waves and P-waves can travel within solids, hence also called body waves. Rayleigh waves generated by the wheel-rail impact have been reproduced and discussed in detail in [50] . Here we show the simulated body waves generated by the wheel-rail impact and their reflections by the sleeper in Fig. 17 . It can be seen that away from the contact point, the velocity is approximately constant across the rail section (see e.g. the section indicated by the dashed red box), meaning the assumption of rigid cross section of the beam model may reasonably apply. However, near the contact point, the beam model is unable to capture the waves propagating from the rail top to bottom, as well as the reflected waves by the sleeper.

Effects on M1 magnitude
Whether the wave propagation should be considered depends on the problem in concern. For example, Rayleigh waves influences the contact solutions as they are generated by the wheel-rail creepage within the contact patch and subsequently propagate through the contact patch [50] . More relevant to this paper are the body waves. For the continuum model, the effective inertia of the rail that participates in the vibration of wheel-track system (e.g. the m eq shown in Fig. 14 ) comes first from the point of contact and then "gradually " expand as the waves spread out, see Fig. 17 . While for the beam model, any vibration always involves the whole cross section, which has a larger inertia. A smaller m eq leads to a smaller M1 peak predicted by the continuum model. It should be noted that although the Winkler bedding model yields nearly identical M1 magnitude to the continuum model, it tends to underestimate the peak force compared to the Kalker's variational method [ 47 , 52 ]. This means there should be other factors that cause the discrepancies of the M1 magnitude between the two models, such as the effect of wave propagation discussed above. In general, taking into account stress wave propagations will lead to a smaller M1 peak. However, the quantitative effect of the wave propagation on the M1 magnitude depends on many factors and thus needs further investigations.

Effects on f 3 resonance
The wave propagation is asymmetrical to the two sides of the contact point; see the asymmetrical velocity field in Fig. 17 . The stress waves travel more freely in the direction away from the support while are more decayed in the direction towards the support. This asymmetry mainly concerns the frequency region around the pin-pin resonance [66] . This may explain why the change of f 3 is asymmetrical about the sleeper support in Fig. 15 .
For the continuum model, the body waves propagate more freely near the rail surface while are more easily reflected and attenuated by the fastenings at the rail bottom (see Fig. 17 ). In contrast, the wave propagation in the beam model does not distinguish between the rail top and bottom. As a result, the beam model might experience more reflected waves at the wheel-rail contact, which further leads to larger frequency fluctuations at f 3 , especially when the contact is near the sleeper support (see Fig. 15 ).

Effects on damping
The damping of the track at higher frequencies (above 1000 Hz) is mainly controlled by the railpad damping. In both the continuum and beam model presented in this paper, the railpads are viscously damped. This means the amount of damping depends on the velocities of the rail and sleeper at rail seats. Fig. 18 compares the rail and sleeper velocities at the rail seat at 6.6 m. It can be seen that the rail velocity magnitude calculated by the beam model is closer to that by the continuum model at the rail top, while is larger than that at the rail bottom. For the sleeper velocity, the beam model also predicts a larger magnitude than the continuum model. The major frequency of the velocities for both models is at around 1100 Hz, which is higher than the pin-pin resonance. Consequently, the continuum model show a lower damping than the beam model at f 3 .

Conclusions
We compare the simulation results of a continuum and a beam FE model for the wheel-rail impact force at a typical rail squat defect. We also compare the simulations with field measurement and observations, which suggests the continuum model is more accurate than the beam model.
The impact force consists of a forced vibration peak M1 followed by free vibration related oscillations with three dominant frequencies f 1 (340 Hz), f 2 (890 Hz) and f 3 (1120 Hz). The three frequencies are independent of wheel speed and defect geometry. They correspond to the eigenfrequencies of the wheel-track system according to the proposed frequency domain model.
The beam model with a Hertzian contact spring overestimates the M1 peak of the impact force. The discrepancy can be reduced but not entirely eliminated by using the Winkler bedding model, because it can better model the wheel-rail contact solution in the longitudinal direction.
Different from the conventional belief that the Timoshenko beam is accurate up to about 1500 Hz for the rail model, we show that the beam model is only accurate up to about 800 Hz in terms of the track receptance. As to the impact force, the valid frequency range of the beam model is further reduced. The beam models produce larger FFT magnitudes at the first characteristic frequency at about 340 Hz. The beam models are only comparable to the continuum model for long and shallow squats, e.g., longer than 50mm while not deeper than 0.3 mm in the case considered in this paper.
The two track models show different damping behaviour around the pin-pin resonance, i.e., between about 800 Hz and 1800 Hz. The damping of the continuum model is larger below and at the pin-pin frequency, whereas the damping of the beam model is larger above the pin-pin frequency. The differences of the damping are caused by the different modelling of the stress wave propagation in the rail. As a result, the contact force dominates at f 3 (1120 Hz) for the continuum model while at f 2 (890 Hz) for the beam model.
The propagation of the body waves in the rail caused by the wheelrail impact is reproduced by the continuum model. We show that the stress wave propagation contributes to the smaller M1 peak and smaller track damping after the pin-pin resonance in the continuum model.
The findings contribute to a better understanding of the dynamic characteristics of the wheel-rail impact force, as well as the root causes for the different simulation results between the two models. In engineering practice, this study can assist engineers in choosing the appropriate assumptions for the wheel, contact and track models when solving wheel-rail impact problems.

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.