Simulation of vertical dynamic vehicle–track interaction using a two-dimensional slab track model

ABSTRACT The vertical dynamic interaction between a railway vehicle and a slab track is simulated in the time domain using an extended state-space vector approach in combination with a complex-valued modal superposition technique for the linear, time-invariant and two-dimensional track model. Wheel–rail contact forces, bending moments in the concrete panel and load distributions on the supporting foundation are evaluated. Two generic slab track models including one or two layers of concrete slabs are presented. The upper layer containing the discrete slab panels is described by decoupled beams of finite length, while the lower layer is a continuous beam. Both the rail and concrete layers are modelled using Rayleigh–Timoshenko beam theory. Rail receptances for the two slab track models are compared with the receptance of a traditional ballasted track. The described procedure is demonstrated by two application examples involving: (i) the periodic response due to the rail seat passing frequency as influenced by the vehicle speed and a foundation stiffness gradient and (ii) the transient response due to a local rail irregularity (dipped welded joint).


Introduction
The use of slab track structures for high-speed railway lines has increased over the last decades [1,2]. By using a slab track design, the need for maintenance can be reduced (compared to ballasted track), which is essential due to the trend of growing traffic and reduced time slots for track works. Other advantages with slab tracks are increased service life, high lateral track resistance and eliminated problems with degradation of ballast. There are, however, also disadvantages with slab tracks, which include higher construction cost, a limited possibility of readjustment of the rail if there is an irregularity in track geometry due to subgrade settlement, and lower vibration/noise absorption [3].
Accurate simulations based on a comprehensive model of the dynamic vehicle-track interaction are required for the design of railway track. Today, there is a large span of existing models that can be used depending on the purpose of the simulation, cf. Connolly et al. [4] and Knothe and Grassie [5]. One of the first models was presented in 1982 by Grassie et al. [6]. In this relatively simple model, a point mass moving on a continuously supported and infinitely long beam is considered, and linearised Hertzian theory is used for the wheel-rail contact.
Imperfections in the vehicle and track may cause severe damage, independently of the track design. Typical examples of imperfections are wheel flats, rail corrugation and differential ballast/soil settlement. In order to understand how these imperfections affect the deterioration of the vehicle and track, simulation of coupled vehicle-track dynamics is needed. Combining the dynamics of the vehicle and track in the modelling procedure started in the early/mid-1990s [7][8][9], albeit for ballasted tracks. A common approach when considering coupled vehicle-track dynamics is to model the track using finite elements together with a modal analysis. Lin and Tretheway [10] applied this approach to model a mass-spring-damper system moving on an elastic beam. This framework was extended with a complex-valued modal analysis technique by Nielsen and Igeland [7].
In the last decade, the research regarding slab track structures has grown. For a floatingslab track system, Li and Wu [11] investigated how load transmission to the soil depends on the length of the slabs. Bezin et al. [12] used a multi-body system analysis to simulate key differences between slab track and ballasted track. Zhai et al. [13] developed a three-dimensional model to investigate the overall vehicle-track system, both for ballasted track and for slab track. Their work included numerical implementations and experimental validation. The concrete slabs were modelled as elastic rectangular plates, and asymmetrical vertical irregularities on the two rails were considered. Lei and Wang [14] modelled the vehicle and track with finite elements and used a moving reference frame (the track model was assumed to be invariant along the track structure). In doing so, the vehicle acts at the same position on the rail throughout the simulation, and the study of the dynamic interaction between a vehicle and a continuous slab track can, therefore, be done efficiently. For discontinuous slabs, Zhang et al. [15] calculated the dynamic wheel-rail contact force due to a sinusoidal rail imperfection. In this work, linear Hertzian theory for the wheel-rail contact was considered, and all simulations were performed in the frequency domain. In recent studies, Sadeghi et al. [16,17] extended the model by including nonlinear properties of the wheel-rail contact and simulations in three dimensions.
In this paper, the vehicle-track interaction is modelled with a complex-valued modal superposition technique for the track and an extended state-space vector approach. The vertical dynamic response can then be calculated by considering a generic initial-value problem. This model was initially developed by Nielsen and Igeland [7] (although for a ballasted track). The performance of the model is demonstrated by considering imperfections in both the track and soil. The track imperfection is a dipped rail weld, whereas the soil imperfection consists of a linear gradient in foundation stiffness and damping.
The paper is organised as follows: Section 2 gives a review of existing slab track systems with reference to the European standard. Section 3 describes the track and vehicle models that are used in this paper as well as the generic initial value problem and its solution procedure. Section 4 presents numerical examples where the rail receptance for different track forms are compared, and the track response due to periodic and transient excitation mechanisms are studied. Finally, Section 5 concludes the paper with a summary and an outlook to future work.

Slab track systems with reference to the European standard
Existing slab track systems can be divided into two major categories, cf. [18], which are continuous rail support systems and direct rail support systems. Continuous rail support systems, where the rail is embedded in an elastomeric layer, were developed from tram applications [2]. When using continuous rail support systems, a settlement-free soil is required due to a minimal possibility of rail readjustment. Direct rail support systems, where the rail is supported at discrete, repetitive locations, is the most common category for high-speed applications and, therefore, the modelling of continuous rail support systems is not included in this paper.
Direct rail support systems can be divided into systems with or without sleepers [18]. When sleepers are used in slab track systems, they can either be laid on the slab, e.g. STEDEF, or embedded in the slab, e.g. Rheda 2000. When sleepers are not used, either prefabricated slabs, e.g. Shinkansen, ÖBB-Porr and Max Bögl, or monolithic systems (paved-in track on civil structures) are used; most of these systems are described in [18]. Independently whether sleepers are used or not, a strengthening layer is placed underneath the slab. If the track is installed in a cold region, a frost protecting layer is required between the soil and the strengthening layer.
The general and specific requirements for different types of slab track systems are described in the European standard prEN16432 [19]. 1 In particular, it is stated that the rail supporting structure (prefabricated elements and/or pavement) shall be designed to distribute the load to the substructure by bending. The calculation models for the load distribution along the track given in [19] assume continuous properties of the slab structure resting on a Winkler foundation with a given bed modulus. Hence, the interaction of the bound and unbound layers is simplified to a model of a track structure with an equivalent thickness that is supported by a representative bed modulus. Two principle designs of the track structure, given in [19], are modelled according to Figure 1(a,b).
There is, however, no indication in [19] regarding which parts of the track that should be designed to distribute the load by bending when prefabricated elements are used. In particular, it is stated that the reinforcement bars in the roadbed are mainly installed to handle temperature variations and not to distribute the vehicle load by bending. Further, requirements for the materials used in the substructure are given in Chapter 10.4 in [19]. Since such components have limited tensile strength and thereby cannot distribute loads  Figure 1. Profile of two typical designs used for slab track systems described in [19]. Track systems correspond to (a) German design with hydraulically bound layer and (b) Austrian design with load distribution slab. by bending, it is evident that they are intended to act only as members that transfer vertical loads in the track structure. Thus, it is not clear which of the track components that should handle the design principle of distributing the load by bending.
In this paper, both the slab track systems sketched in Figure 1 will be modelled. Depending on the physical properties of the different layers, different models can be used. Throughout this paper, the rail and concrete panels are modelled using finite elements, while the foundation/soil and the frost protection layer are combined into an effective stiffness and damping of a Winkler foundation. For the slab track type in Figure 1(a), the roadbed is also included in the Winkler foundation, while the roadbed is modelled with finite elements for the slab track type in Figure 1(b).

Modelling and simulation of vertical vehicle-track interaction
Two different types of slab track models are considered, see Figures 2 and 3. In Figure 2, the slab is modelled by one continuous layer of beam elements, whereas in Figure 3, the slab is modelled by two layers of beam elements. In the two-layer model, the upper layer containing the discrete slab panels is described by (possibly coupled) beams of a given length p. In both models, the rail is discretely supported at rail seat distance L = 0.65 m. The bottom beam layer (panel for the track model shown in Figure 2 and roadbed for the track model shown in Figure 3) is supported by non-interacting springs and dampers (Winkler foundation) with the bed modulus k f (X) = 100 MN/m 3 and viscous damping c f (X) = 82 kNs/m 3 . The corresponding foundation stiffness per unit beam length is obtained by multiplying the bed modulus with the width of the slab. For the track model shown in Figure 3, the length p of the discrete panels is 5.2 m (which corresponds to eight rail seat distances) and the discrete connection between two panels is centred between two adjacent rail seats.
The slab track model in Figure 2 represents the type of slab track illustrated in Figure 1(a). In this model, it is assumed that the layers beneath the continuous concrete panel do not distribute any loads by bending. Hence, it is sufficient to model the rail and concrete panel as two layers of beam elements, while remaining layers are incorporated  . Sketch of track model where the slab is modelled as two layers of beam elements. The model contains three layers of beams: rail (r), discrete panels of concrete slab (s 1 ) and continuous concrete roadbed (s 2 ). The concrete base is supported by a Winkler foundation. Note that the coupling between layers (s 1 ) and (s 2 ) is continuously distributed. This slab track model is used to model the slab track type in Figure 1 in the Winkler foundation. Similarly, the slab track model in Figure 3 represents the slab track illustrated in Figure 1(b). Here, the reinforced roadbed may distribute the load by bending and, therefore, the rail, concrete panels and reinforced roadbed are all modelled as layers of beam elements, while the frost protection layer and soil are incorporated in the Winkler foundation. Note that the spring-damper coupling between the two concrete layers is continuously distributed.
The load from the vehicle is assumed to be symmetrically distributed between the two rails and, therefore, only half of the slab and one rail need to be considered. Moreover, only vibrations in the vertical direction are studied in this paper.
In the simulations, the moving wheel-rail contact forces are converted into spatially stationary nodal forces on the finite element track model. The conversion implies that the true local deflections under the contact forces are somewhat underestimated between nodes [20]. This effect can be reduced by increasing the number of beam elements per unit length. In the current study, the rail and concrete slabs between two rail seats are each modelled by 16 elements. For such a discretisation, the effect is negligible. Convergence studies in terms of the influence of the number of beam elements on receptance, wheel-rail contact forces and bending moment in the panels have been performed. From these studies, it was found that convergence of the track responses could be seen already when four elements were used. 2 However, 16 elements per rail seat distance are preferred to improve the load representation without introducing unmanageable computational times.

Railway track model
In this study, finite elements are combined with a complex-valued modal analysis to determine the dynamic response of the track. The slab and rail are both modelled using the Rayleigh-Timoshenko beam elements. The disadvantage of using finite elements in simulations of the vehicle-track dynamics (compared to multi-body simulations) is that the system of equations of motion becomes large and coupled. The computational effort can, however, be reduced by using a modal superposition technique that decouples the equations of motion.
Throughout this paper, the same rail properties are used and hence the rail is labelled with index 'r' for both track models. In Figure 2, the panel is labelled 's', whereas in Figure 3, the panel and roadbed are labelled 's 1 ' and 's 2 ', respectively. Each layer of beam elements has bending stiffness EI A (X), 3 shear stiffness kGA A (X), mass m A (X) per unit beam length, rotatory inertia m A r 2 A (X) per unit beam length, width b A and height h A (X) that may vary along the longitudinal track coordinate X. However, in the present paper all beam properties of each layer are taken as constant. The connection between each vertical pair of adjacent nodes in the different layers is modelled as a spring and viscous damper in parallel and denoted as B C/D . 4 For example, the rail pad stiffness for the track model in Figure 2 is denoted k r/s . For the rail pad connection, also rotational stiffness k r/s δ 2 /12 is accounted for, where δ = 0.15 m is the width of the pad. The values used for the track in the numerical simulations are presented in Tables 1 and 2.
As indicated in Figures 2 and 3, the slab and rail are clamped at both ends. Since the model uses finite elements and the analysis is in the time domain, a finite length of track is suitable. This implies that reflections of vibrations will occur at the ends. For this reason, structure-borne vibration energy is not fully transmitted away from the excitation. This effect will, however, be negligible at the centre of the track if the length of the model is sufficient. Here, the length of the track model is set to 78 m (corresponding to 120 rail seat distances), which results in a small boundary effect on the track receptance, see Section 4.1. However, in the simulations of vehicle-track interaction, these boundary effects are small and considered negligible. 5 Since the rail is modelled by Rayleigh-Timoshenko beam elements, there is a (nonphysical) discontinuity in slope over the element boundaries due to the applied interpolation polynomials, cf. [21]. In the numerical examples, this can be observed to cause small Table 1. Beams properties in the track models.
Concrete panel, Figure 2 Concrete panel, Figure 3 Concrete roadbed, Figure 3 Rail, Figures 2 and 3 Table 2. Stiffness and damping of resilient layers in the track models.
Stiffness Damping disturbances in the wheel-rail contact force each time an element boundary is passed (which would not be present if Euler-Bernoulli beam elements were used). In the numerical example in Section 4.2, this effect is rather pronounced since no wheel/rail irregularity is considered. Mass, stiffness and damping matrices for the track structure are derived using thirddegree polynomial shape functions derived for Rayleigh-Timoshenko beam theory, cf. [21,22]. Let x t (t) be a vector containing the degrees-of-freedom (DOFs) of the track model and let F t (t) contain the external loads (in particular, the wheel-rail contact force F a (t)). The coupled equations of motion can then be written in state-space form (superscript t denoting track) as where where K t , C t and M t are the stiffness, damping and mass matrices of the track, respectively. By utilising that A t and B t are symmetric, which makes the left and right eigenvectors each other's transpose, the complete modal solution to Equation (1a) is given by where ω n are angular eigenfrequencies, ρ (n) are eigenvectors and I denotes the unit matrix. 6 By solving the eigenvalue problem expressed in Equation (2), the modal matrix, P, can be assembled. The modal matrix works as a mapping between the spatial and modal domains as where q t (t) is the modal displacement and Q t (t) is the modal load vector. By using the orthogonality of the modal matrix, the decoupled equations of motion can be written as where diag(a n ) When q t (t) and Q t (t) have been calculated (see Sections 3.3 and 3.4), it is straightforward to determine the physical displacements, velocities and accelerations. For more general information about the track modelling approach, see [7].

Vehicle model
As indicated in Figures 2 and 3, the vehicle model consists of one bogie. More advanced vehicle models, e.g. including several bogies and car bodies, can be implemented. Here, a single bogie model is considered since the secondary suspension acts as a dynamic filter isolating the car body from the bogie in the frequency range where the dynamics of the track is significant (20 < f < 1500 Hz) [5]. 7 The bogie is modelled by six DOFs and is containing wheel-rail contact stiffness and primary suspension stiffness k and damping c. Four of the vehicle DOFs represent the motion of the vehicle and bogie, while the two remaining DOFs are interfacing the track and are used in the formulation of the constraint equations.
The weight of half of the car body is accounted for by a static point-force acting at the centroid of the bogie. The axle load of the vehicle is 17 tonnes. The bogie also contains two unsprung masses M b1 and M b2 , corresponding to the two wheelsets, separated by distance from each other, and a bogie frame with mass M b3 and inertia J b3 . The parameter values for the vehicle are presented in Table 3 and were collected from [23] (except M 0 ).
The DOFs of the vehicle are collected in two vectors: x v a = {x a1 x a2 } T that contains the two massless DOFs that are interfacing the track and x b3 x b4 } T that contains the non-interfacial DOFs (representing the two wheelsets and one bogie). Superscript v denotes vehicle and will be used throughout this paper. The equations of motion for the vehicle can then be written as where F a (t) are contact forces between the wheels and the rail and F ext b contains all external forces (in this study only gravity loads). The stiffness matrix K v includes the Hertzian contact stiffness, which is given by where C H is the Hertzian constant and the Macaulay brackets are defined as

Coupling of vehicle and track models
In order to formulate constraint equations, consider the four third-degree interpolation polynomials, 8 cf. [22], where β j = 12EI/(kGAl 2 j ) and ξ j ∈ [0 l j ] is the local coordinate of element j with length l j in the finite element model of the rail. The constraint between the vehicle and the track can then be formulated as [24] x where x irr contains the prescribed wheel/rail surface irregularity, P int is the time-variant partition of the modal matrix corresponding to the rail DOFs that are adjacent to the current position of the vehicle DOFs and N is a block matrix containing the interpolation polynomials given in Equation (7). The constraints on the interfacial velocities can then be calculated as the time derivative of Equation (8) asẋ Similarly, the constraints of the interfacial accelerations are obtained by a second differentiation asẍ v where Here, R contains the Coriolis and vertical accelerations and S includes the effect of the centripetal acceleration and the longitudinal acceleration of the vehicle. Moreover, the modal load vector is obtained as Further details on the derivation of the constraint equations are given in [24].

Solution of the vehicle-track interaction problem
In order to solve the interaction problem, the mixed extended state-space vector z is defined as whereF a = F a (t)dt. The total coupled and time-variant system is written as where An initial-value problem can now be formulated aṡ By exploiting that the modal analysis yields complex-conjugated sets of modal parameters, the real and imaginary parts of Equations (4a) and (9)-(12) can be separated making the initial-value problem in Equation (16) real valued. Moreover, it is noted that both A and B are sparse matrices.

Calculation of track response including the influence of residual terms
In the numerical examples (see Section 4), typical outputs of interest are wheel-rail contact forces, bending moments in the slab and load distributions on the foundation. The contact forces are included inż and can be determined from the initial-value problem. However, in order to calculate a track response, additional calculations are required to compensate for the truncated modes. From the solution of the initial-value problem, physical displacements can be obtained by multiplying the modal displacement with the modal matrix. In this paper, the influence of the truncated modes on the total physical displacement is added as a quasistatic contribution, cf. [25]. By assuming that the dynamic effects of the truncated modes can be neglected, Equation (4a) is reduced to where the subscript τ indicates that the equation is only valid for the truncated modes. The contribution from the truncated modes to the physical displacements and velocities (denoted δ τ x t (t) and δ τẋ t (t), respectively) can then be expressed as Exploiting that the flexibility elements δe jk can be calculated from a frequency-dependent flexibility matrix implies that where r is the number of retained modes and e jk (ω ref ) is the flexibility matrix at a much lower frequency than any of the eigenvalues of the truncated modes. The sectional forces (shear force and bending moment) at a given node in the track model are calculated using the corresponding element stiffness matrix, F, and the associated physical displacements and rotations. Let the subset of x t containing the physical displacements and rotations (including the quasistatic contribution from the truncated modes) at the two nodes of the given element be denoted as n. The vectorÑ(t), containing the time-variant sectional forces, is then calculated as Fn(t). (20) The bending moment for the given node is obtained by extracting the associated term iñ N. To calculate the load distribution on the foundation, the physical displacements and velocities of the bottom beam layer (panel for the track model shown in Figure 1(a) and roadbed for the track model shown in Figure 1(b)) are extracted from x t andẋ t , respectively. For a given longitudinal coordinate X, the time-variant load L d (X, t) per unit beam length [N/m] is calculated as where k f (X) and c f (X) are the properties of the foundation at the given element, b is the width of half of the slab (since only half of the slab is used in the track model) and n contains the physical displacements and rotations (including the quasistatic contribution from the truncated modes) at the two nodes adjacent to the coordinate X. The four interpolation polynomials in N are evaluated at the corresponding local coordinate ξ of the corresponding element.

Numerical examples
In this paper, three demonstration examples are presented. In Section 4.1, the track response is investigated by calculating the rail receptance for the two slab track models, and comparing the results with the receptance of a traditional ballasted track model. In the other two demonstration examples, the track model with discrete panels in Figure 3 is used. In Section 4.2, the influence of a foundation stiffness gradient on the wheel-rail contact force, the corresponding stiffness gradient at the rail level and the load distribution on the foundation is investigated. Finally, the effect of a dipped welded rail joint is investigated in Section 4.3. The weld is modelled by prescribing a vertical rail irregularity (accounted for in x irr , see Equation (8)). The resulting wheel-rail contact force and bending moment in the slab are evaluated for different depths of the dip. Furthermore, the influence of the weld location, relative to the ends of the discrete panels, is investigated. The demonstration examples in Sections 4.2 and 4.3 show significant dynamic effects. The magnitudes of these dynamic effects, which would not have been captured by calculations according to the European Standard [19] (which only applies static calculations), are easily analysed using the presented model.

Track receptance
For a comparison of the different track forms, receptances of the track models are calculated in the frequency domain. The harmonic force excitation is applied on the rail and the response is calculated on both the rail and the panel. Let x t r denote a (complex-valued) displacement on the rail, x t s be a displacement on the concrete panel and F t r be a force applied on the rail. The considered receptances are then given by In the presented numerical example, the excitation and the response are evaluated at a rail seat. Magnitudes and phases of the complex-valued receptances H t r and H t s are presented in Figures 4 and 5.
For the ballasted track model, the same rail and soil properties are used, but the slabs/panels are replaced by discrete sleepers. Each sleeper is modelled by Rayleigh-Timoshenko beam elements. The pad stiffness and viscous damping are higher for the ballasted track (k p = 120 kN/mm, c p = 25 Ns/mm), compared to the slab track models since this is common design practice.
From Figure 4, it can be seen that all track types have two resonances and two antiresonances in common (in the studied frequency range). At the first resonance frequency (40-140 Hz depending on track model), the slab/sleeper and rail move in phase, while they move out of phase at the second resonance frequency (160-300 Hz) [5]. The anti-resonance at about 950 Hz is the so-called pinned-pinned frequency, where the rail vibrates with a wavelength corresponding to the distance between two rail seats while the slab/sleeper does not move. The small difference in magnitude and frequency of the 'pinned-pinned' mode between the slab track models and the ballasted track model is due to differences in the vertical and rotational rail pad stiffness and damping. The point receptance on the rail between two rail seats is similar to the response on the rail at the rail seat up to about 500 Hz (not shown here). Above 500 Hz, the receptances differ, and at the 'pinned-pinned' frequency, there is a resonance peak. The presented receptances can be compared with the calculated receptance for a floating-slab track (up to 500 Hz) [11]. As expected, the same resonances can be identified, but with slightly different magnitudes and frequencies due to different input data. In Figure 5, H t s is presented for the different track models. The receptance varies significantly between the different models. However, it can be seen that the first resonance observed in H t r (between 40 and 140 Hz) has a pronounced effect also on H t s . Moreover, the resonance at about 850 Hz for the continuous one-layer model can be seen as a small peak in Figure 4(a). It is important to emphasise that the presented results are receptances based on the input values presented in Tables 1 and 2. Although the input values are selected to be as realistic as possible, no experimental calibration/validation has yet been performed for the slab track model. The ballasted track model has been verified and compared with physical testing, cf. [26].

Foundation stiffness gradient
For the track model in Figure 3, the influence of a foundation stiffness gradient (per unit beam length) on the wheel-rail contact force and the dynamic response of the track is investigated. In this case, the dynamic excitation is dominated by the rail seat passing frequency (periodic excitation) and a transient excitation due to the stiffness gradient. The foundation stiffness is altered along the track by a stepwise change in stiffness and damping of the Winkler foundation in adjacent beam elements according to Figure 6. Based on the stiffness matrix of the track model, the static stiffness at the rail level along the track is calculated for the given foundation stiffness gradients, see Figure 7. As expected, the stiffness at the rail level has a periodic pattern due to the discrete rail supports. For the given slab track design, it can be evaluated that gradients in the foundation stiffness of 0, 5 and 15 (kN/mm)/m 2 correspond to average stiffness gradients of 0, 0.5 and 1.5 (kN/mm)/m at the rail level. 9 Figures 8-10 show the calculated wheel-rail contact forces for three different vehicle speeds. For a given foundation stiffness gradient, it is observed that the influence of the vehicle speed on the amplitude of the wheel-rail contact force due to the periodic excitation by the rail seat passing frequency is significant. For v = 100 km/h, the two positive stiffness gradients lead to an increase in contact force amplitude, see Figure 8, whereas the same gradients lead to a decrease in contact force amplitude at v = 150 km/h, see Figure 9. For v = 300 km/h, the influence of the studied stiffness gradients is negligible, see Figure 10.
In the results illustrating the calculated wheel-rail contact forces, a low-pass filter has been applied for aesthetic reasons since Rayleigh-Timoshenko beam theory has been used when deriving the finite element matrices. As discussed in Section 3.1, in Rayleigh-Timoshenko beam theory, there is a (non-physical) discontinuity in slope over the element boundaries due to the applied interpolation polynomials, cf. [21]. Without a low-pass filter, this can be observed to cause small disturbances in the wheel-rail contact force each time an element boundary is passed. The cut-off frequency is selected to be below the element passing frequency, which for the lowest considered speed (100 km/h) is 684 Hz. Therefore, the cut-off frequency was set to 650 Hz. It has been confirmed that the influence of higher frequencies on the dynamic response (except the spurious effect of the element boundaries) is negligible.
In Figures 8-10, it is observed that the wheel-rail contact force is dominated by the periodic response due to the rail seat passing frequency. For direct rail support systems, the rail seats induce a parametric excitation due to the repetitive variation in track stiffness. Therefore, the track dynamics at the rail seat passing frequency will, to a large extent, determine the magnitude of the wheel-rail contact forces. For stationary harmonic excitation, the amplitude of the contact force is proportional to the inverse of the system receptance, cf. [6,27]. Thus, the influence of the vehicle speed on the wheel-rail contact force amplitude can be explained by the receptance of the vehicle-track system. By letting x v be the displacement of the wheelset/unsprung mass and F v the force applied on the wheelset/unsprung mass, the receptance of the wheelset/unsprung mass is given by (neglecting the influence of the primary suspension and the sprung components of the vehicle) The system receptance is obtained by adding the receptance of the vehicle and track as 10 For the two-layer slab track model, Figure 11 shows the magnitude of the system receptance for three different soil stiffnesses (low, mean and high stiffness when the stiffness gradient is 15 (kN/mm)/m 2 according to Figure 6). The vertical lines indicate the rail seat passing frequency when v = 100 km/h, v = 150 km/h and v = 300 km/h. As can be seen from the figure, at vehicle speed v = 100 km/h, the rail seat passing frequency coincides with a local minimum in the system receptance for the slab track on the high foundation stiffness. Further, at this frequency, the slab track on the high foundation stiffness has the lowest system receptance for the three foundation stiffnesses considered. This explains why the amplitude of the wheel-rail contact force increases with increased foundation stiffness when v = 100 km/h, see Figure 8. An opposite situation occurs at v = 150 km/h, where the track model on the low foundation stiffness leads to the lowest system receptance, which explains the negative slope in the amplitude of the wheel-rail contact force in Figure 9. Moreover, it can be seen that the system receptance is independent of the stiffness of the soil at v = 300 km/h, which explains why also the amplitude of the wheel-rail contact force is not influenced by the stiffness gradient in Figure 10. Figure 12 shows the system receptance, but for the ballasted track model. From the figure, it can be seen that there is a difference in the receptance for the different foundation stiffnesses up to much higher frequencies, compared to the slab track model. This phenomenon is expected since the panel, together with the roadbed, is a much stiffer structure compared to the sleepers. Finally, Figure 13 shows that an increase in foundation stiffness gives an increase in the wheel-rail contact force amplitude when v = 150 km/h for the ballasted track, which is in line with the conclusions drawn from Figure 12. Based on the calculated displacement and velocity fields of the roadbed, the influence of the foundation stiffness gradient on the load distribution at the interface between roadbed and foundation has been investigated. Figure 14 shows the load distribution on the foundation at a time instant for three different cases with uniform soil stiffness (low, mean and high stiffness when the stiffness gradient is 15 (kN/mm)/m 2 according to Figure 6). In all simulations, the thickness of the roadbed is h s 2 = 20 cm and the vehicle speed is v = 300 km/h. For the presented set of track model input data, it can be seen that the soil stiffness has a low impact on the magnitude and extent of the load distribution (less than 1%). Even though the effect is small in the presented figure, a higher foundation stiffness leads to a higher maximum value and a more narrow distribution compared to the case with a lower foundation stiffness. 11 A stiffness increase implies a displacement decrease, which (for the considered vehicle model) leads to similar load distributions independently of soil stiffness. Note that the results in Figure 14 are only considering half of the slab. Note that the shape and peak value of the load distribution depends on the thickness of the roadbed. A thicker roadbed would imply a wider distribution and a reduced peak value. If the roadbed was decreased significantly, e.g. h s 2 = 5 cm, two peaks (one from each wheelset) would have been seen, and the influence of the discrete panels had been more pronounced (not shown in Figure 14). For results illustrating the influence of the thickness of the roadbed, see [28].
For each track position, and for the three different stiffness gradients defined in Figure 6, the maximum value of the load distribution due to the moving vehicle load has been identified and is presented in Figure 15. 12 The vertical lines indicate where one discrete panel ends, and the next panel begins. Two periodic patterns are visible in the figure, one is due to the discrete rail seats, and one is due to the discrete panels. Further, the foundation stiffness gradients start at track position x = 35.9 m and end at track position x = 41.4 m. These track positions correspond to minima and maxima in the load distribution on the foundation. These extreme points occur since the physical displacement field is smooth  over the start and end points of the stiffness gradient, while the foundation stiffness has discontinuous derivatives at these locations, cf. Figure 6 and Equation (21).

Dipped welded rail joint
In a third demonstration example, the influence of a dipped welded rail joint leading to a more transient track response is illustrated. The imperfection caused by the weld is modelled as a purely geometrical imperfection by symmetrically prescribing the vertical level of the rail in the area around the weld. The symmetrical imperfection consists of two second-order polynomials as, cf. [29], where X is the track position,X is the centre position of the weld, d is the depth of the weld and D is the length of the rail that is affected by the irregularity. An illustration of the irregularity with depth d = 1 mm is shown in Figure 16. In all simulations in this section, the vehicle speed is 300 km/h, and the weld is located at X = 38.675 m, which is above the centre of a panel (and therefore also midway between two rail seats). Figure 17 shows the wheel-rail contact force at the leading wheelset as a function of track position. From the figure, it can be seen that the depth of the weld has a large influence on the maximum value of the wheel-rail contact force. For the given input data, the maximum value of the wheel-rail contact force is increased by 46% for d = 0.5 mm and by 93% for d = 1 mm. Apart from the large dynamic loads when the leading wheelset passes the weld atX = 38.675 m, large dynamic loads can also be seen aroundX = 41.175 m, which is the position when the trailing wheelset passes the weld. The impact load from the trailing wheelset has a significant effect also at the leading wheelset due to the low stiffness and damping of the rail pads which imply a low decay rate of rail vibrations.
The effect of a dipped welded rail joint is further studied by calculating the bending moment in the panels (following the procedure described in Section 3.5). The bending moment is generally higher below the rail seats compared to between the rail seats, and, therefore, the bending moment is calculated in the panel below the first rail seat after the weld. In Figure 18, it is observed that the depth of the weld has a large impact also on the magnitude of the bending moment. Note that the results in Figure 18 are only considering half of the slab.
Finally, the influence of the weld location, relative to the end of the discrete panels, is investigated. Several simulations were performed where the centre position of the dipped welded rail joint was changed. For each position of the weld, the bending moment in the panel was calculated below the position of the weld and at the adjacent rail seats. The maximum bending moment of the considered simulation was then identified as an indicator for the considered weld position. Figure 19 shows the maximum bending moment as a function of the position of the weld for three different panel lengths, p. The considered panel lengths are p = 2.6 m, p = 5. figure, it can be seen that there is a periodic pattern due to the discrete panels. In particular, the bending moment is always lower at the joint between two panels compared to other positions (independently of panel length).

Conclusions
A general procedure for analysing the vertical dynamic vehicle-track interaction has been presented and implemented successfully in a slab track context. The performance of the procedure has been illustrated by demonstration examples that contain structural irregularities. Based on the demonstration examples, it was found that the influence of a foundation stiffness gradient on the amplitude of the wheel-rail contact force is vehicle speed dependent due to the rail seat passing frequency and the excitation of the fundamental vehicle-track system resonance. For the given slab track models, it was observed that the influence of the studied foundation stiffness gradients on the wheel-rail contact force was negligible when high speed was considered (v > 200 km/h). Even if the amplitude of the wheel-rail contact force was found independent of the foundation stiffness gradient at high speed in the slab track context, other dynamic responses such as the bending moment in the slab and the load distribution on the foundation were affected. Moreover, it was found that a geometrical rail imperfection, here exemplified by a dipped welded rail joint, has a large impact on the wheel-rail contact forces as well as on the bending moment in the panels. In addition, it was found that the spatial location of the weld, relative to the transition between two slab panels, affects the results significantly. Thus, in the described demonstration examples, significant dynamic effects were observed. In contrast to the European Standard [19], which only applies a static calculation approach, the magnitude of these dynamic effects is analysed with the presented model.
In future work, the model of the slab will be extended to three spatial dimensions, and it will be determined whether it is sufficient to model the panels as Rayleigh-Timoshenko beams or if a three-dimensional (3D) slab track model with plate elements is necessary. Since the computational cost can be controlled by the number of used modes, the presented analysis procedure seems suitable also for a 3D model. Moreover, transition zones between slab track and ballasted track will be investigated to evaluate an optimal design. Notes 1. How to interpret the European Standard is further elaborated by the International Union of Railways [30]. 2. Depending on the accepted error tolerances. 3. Replace A with the label of the considered layer (r, s, s 1 or s 2 ). As an example, the bending stiffness of the rail is denoted EI r (X). 4. Replace B with k or c depending on if it is a spring or a damper and let the spring/damper connect layers C and D. 5. This has been confirmed by a convergence study. 6. Underline indicates complex-valued quantities. 7. When acoustic effects of the track are studied, frequencies up to 5000 Hz is typically studied. 8. The interpolation polynomials are the shape functions derived from Rayleigh-Timoshenko beam theory. 9. Since the foundation is described by a Winkler model where the stiffness is given as stiffness per unit beam length, the unit of the foundation stiffness gradient is [(N/m)/m 2 ]. 10. Here, the influence of wheel-rail contact stiffness on the system receptance is neglected as it is frequency-independent. 11. This effect is more pronounced if only one wheelset is used as vehicle model or if the distance between the wheelsets is increased. 12. Similar to Figure 14, the results in Figure 15 are highly dependent on the thickness of the roadbed.