Transition radiation in an infinite one-dimensional structure interacting with a moving oscillator—the Green’s function method

Abstract Transition zones in railway tracks are areas with considerable variation of track properties (e.g., foundation stiffness) which may cause strong amplification of the response, leading to rapid degradation of the track geometry. Two possible indicators of degradation in the supporting structure are identified, namely the wheel-rail contact force and the power/energy input by the vehicle. This paper analyses the influence of accounting for the interaction between the vehicle and the supporting structure on the contact force and on the power/energy input. To this end, a one-dimensional model is formulated, consisting of an infinite Euler-Bernoulli beam resting on a locally inhomogeneous Kelvin foundation, interacting with a moving loaded oscillator that has a nonlinear Hertzian spring. The solution is obtained using the Green’s-function method. To obtain the Green’s function of the inhomogeneous and infinite beam-foundation sub-system, the finite difference method is used for the spatial discretization and non-reflective boundary conditions are applied. Accounting for the interaction between the moving oscillator and the supporting structure generally leads to stronger wave radiation, caused by the variation of the vertical momentum of the moving mass. Results show that for relatively high velocity and small transition length, the maximum contact force as well as the energy input exhibit a significant increase compared to the moving constant load case. Furthermore, for relatively high velocities, the maximum contact force also increases significantly with increasing stiffness dissimilarity, findings which supplement the existing literature. Finally, the two degradation indicators can be used in the preliminary stages of design to assess the performance of railway track transition zones.


a b s t r a c t
Transition zones in railway tracks are areas with considerable variation of track properties (e.g., foundation stiffness) which may cause strong amplification of the response, leading to rapid degradation of the track geometry. Two possible indicators of degradation in the supporting structure are identified, namely the wheel-rail contact force and the power/energy input by the vehicle. This paper analyses the influence of accounting for the interaction between the vehicle and the supporting structure on the contact force and on the power/energy input. To this end, a one-dimensional model is formulated, consisting of an infinite Euler-Bernoulli beam resting on a locally inhomogeneous Kelvin foundation, interacting with a moving loaded oscillator that has a nonlinear Hertzian spring. The solution is obtained using the Green's-function method. To obtain the Green's function of the inhomogeneous and infinite beam-foundation sub-system, the finite difference method is used for the spatial discretization and non-reflective boundary conditions are applied. Accounting for the interaction between the moving oscillator and the supporting structure generally leads to stronger wave radiation, caused by the variation of the vertical momentum of the moving mass. Results show that for relatively high velocity and small transition length, the maximum contact force as well as the energy input exhibit a significant increase compared to the moving constant load case. Furthermore, for relatively high velocities, the maximum contact force also increases significantly with increasing stiffness dissimilarity, findings which supplement the existing literature. Finally, the two degradation indicators can be used in the preliminary stages of design to assess the performance of railway track transition zones.

Introduction
Transition zones in railway tracks are areas with substantial variation of track properties (e.g., foundation stiffness) encountered near rigid structures such as bridges, tunnels or culverts. These zones require 3-8 times more frequent are identified. The current paper has some similarities to the works of Lei and Mao [26] and of Ang and Dai [19] , but is distinct in the following ways. Firstly, the vehicle models in both studies are more complex than the one used in the current paper, making the observations difficult to relate/compare to the case of a single moving constant load. Secondly, both works are limited to piecewise-homogeneous foundations, while the current paper addresses a smooth stiffness and damping variation of the supporting structure. Thirdly, the infinite extent of the system is properly accounted for in the present work. Finally, the extent of the parametric study presented in both previous works is considerably enlarged in the current study, and thus, more general observations of the system's behaviour can be made, relevant for engineering practice.
Modelling the railway track as a 1-D system (with frequency-independent parameters) has its shortcomings. For example, the wave propagation in the vertical direction as well as surface waves are not captured by this model, aspects which can influence the dynamic response. However, the qualitative observations and the phenomena identified with this simplified model will be present in a more comprehensive model. Therefore, for the purpose of studying the influence of the vehicle's own degree-of-freedom on the transition radiation phenomenon, the 1-D model is sufficient. It must be emphasized that this paper focuses on transition radiation as the process of energy emission, which is characterized by the energy flux though a boundary that encompasses the emitter, with no emphasis on the type of waves that are generated (propagating and/or evanescent).

Problem statement
In this section, a 1-D model is formulated, consisting of an infinite Euler-Bernoulli beam resting on a smoothly inhomogeneous Kelvin foundation, interacting with a moving one-mass loaded oscillator ( Fig. 1 ). The infinite system is divided into three domains: the computational domain (i.e., x ∈ [0 , L ] as can be seen in Fig. 1 ) that contains the entire inhomogeneity of the supporting structure and its vicinity, and two homogeneous semi-infinite domains, one at each boundary of the computational domain. The initial position (i.e., at t = 0 ) of the moving oscillator is at x = 0 and its oscillations are restricted to the computational domain. The equation of motion for the three domains reads x ≥ L, x ≥ L, represent generalized derivatives with respect to space and time, respectively. Compared to the classical derivatives, the generalized ones can have discontinuities, and are necessary in Eq. (1) due to the Dirac delta function on the right-hand side, which is a generalized function [49] . Furthermore, all parameters in Eq. (1) have been scaled by the beam's bending stiffness EI; ρ represents the scaled mass per unit length of the beam; k d , l and c d , l are the scaled (homogeneous) foundation stiffness and damping of the left semi-infinite domain, respectively; k d , r and c d , r represent the same quantities of the right semi-infinite domain; k d , c (x ) and c d , c (x ) are the scaled foundation stiffness and foundation damping of the computational domain, respectively; Q c (t) ,Q 0 and v represent the time-dependent contact force in the computational domain, the steady-state contact force (due to the dead load of the car body) and the velocity of the moving oscillator, respectively; δ( . . . ) denotes the Dirac delta function. Finally, w l ,w r and w c are the unknown displacements of the left and right semi-infinite domains, and of the computational domain, respectively. The space and time dependency of the unknown displacements is omitted from most expressions for brevity. It must be emphasized that this model does not incorporate the influence of the discrete nature of the rail support. Moreover, by adopting ≤ and ≥ signs in the definition of w (x, t ) ,Q (t) ,k d (x ) , and c d (x ) for all domains it is emphasized that there is continuity in these quantities at the interfaces between the three domains.
In this formulation, the vehicle is modelled as a loaded one-mass oscillator (i.e., loaded mass-spring system), where the constant dead load Q 0 (which includes the dead weight of the wheel) corresponds to the vehicle weight on the wheel, the mass M represents the wheel of the vehicle and the spring represents the contact elasticity between the wheel and the rail. The car body and the bogies could be modelled by additional suspended masses. However, the frequencies at which these suspended masses oscillate are much lower than the oscillation frequency of the wheel in the interaction with the rail. Therefore, the degrees of freedom of the suspended masses of the car body and bogies are neglected. In the remainder of the paper, the term oscillator is used to denote the loaded one-mass oscillator. The equation of motion for the oscillator reads where u is the displacement of the mass and overdots denote classical derivatives with respect to t. The interaction between the oscillator and the beam-foundation sub-system is described by the Hertzian contact law, where the possibility of contact loss is incorporated as presented in the following equation: where C H is Hertz's constant and H( . . . ) denotes the Heaviside function. It must be noted that for large transition lengths, the inertia of the suspended masses may have a non-negligible effect on the behaviour due to the low frequencies of the interaction [32] .
From Eq. (1) it can be deduced that the shear force is continuous everywhere except under the moving load where the force is discontinuous. This discontinuity is accounted for in the equation of motion ( Eq. (1) ), for all time moments, due to the presence of the generalized derivatives. Consequently, as interface conditions between the domains, continuity in displacement and slope, as well as in shear force and bending moment is imposed. The set of boundary conditions is completed by imposing that, due to the presence of damping, the displacements of the left and right domains are zero at infinite distance from the moving oscillator: where primes denote classical derivatives with respect to x . The computational domain could be made very large such that the influence of the initial conditions on the response dies out before the oscillator reaches the transition zone. However, for computational efficiency, x = 0 should be as close as possible to the inhomogeneity (i.e., the computational domain should be as small as possible), meaning that the choice of initial conditions affects the response in the transition zone. Prior to reaching a transition zone, the response caused by a train can generally be considered to be in the steady state. Consequently, the initial conditions must be chosen such that the system reaches the steady-state regime instantaneously at the start of the simulation, and are thus based on the steady-state field of the approaching load, also referred to as the eigenfield w e . For the considered problem, in the steady state, the moving oscillator is equivalent to a moving constant load where the magnitude of the force is equal to Q 0 , and the eigenfield for such a system reads [50,51] w e (x, t ) = where A 1 ,A 2 ,B 1 and B 2 represent complex-valued amplitudes and k e 1 ,k e 2 ,k e 3 and k e 4 represent the complex wavenumbers of the eigenfield. Their expressions are given in Appendix A . It must be noted that w e in Eq. (9) is real-valued.
The initial conditions based on the eigenfield are correct only if the initial fields (i.e., displacement and velocity fields) are not disturbed by the inhomogeneity because the eigenfield relates to a homogeneous Kelvin foundation. Consequently, x = 0 should be positioned such that the initial fields based on the eigenfield are negligibly small at the inhomogeneity (they cannot be precisely zero due to the eigenfield's exponential decay in space). The initial conditions of the left domain are imposed as the eigenfield part behind the oscillator, the ones of the computational domain are imposed as the eigenfield part in front of the oscillator, and the initial conditions of the right domain are trivial: The initial conditions of the oscillator should be chosen such that it also reaches the steady state instantaneously at the start of the simulation. Considering that the contact force is constant and equal to Q 0 in the steady state, the initial conditions for the oscillator are given as follows: Eqs. (1) to (13) constitute a complete description of the problem. Next, the solution method is presented.

Solution method
The system described in Section 2.1 is nonlinear due to the nonlinear contact relation ( Eq. (3) ). However, the beamfoundation sub-system behaves linearly. This allows to obtain the solution using the time-domain Green's-function method [48] which is explained next.

General procedure
The goal of the Green's-function method is to express the displacement of the beam at the contact point in terms of the unknown contact force Q (t) . The obtained displacement can be substituted in the contact relation ( Eq. (3) ). Then, the displacement of the mass can be expressed in terms of the contact force, and the resulting expression can be substituted in the equation of motion of the mass. The resulting integro-differential equation can then be solved for the unknown contact force numerically. This approach is suited for any type of vehicle model. However, taking advantage of the linearity of the vehicle model used in this paper (except for the contact spring), the displacement of the vehicle is expressed in terms of the unknown contact force through the Green's-function method too. Then, the contact relation reduces to an integral equation for the contact force, and is solved iteratively.
To express the displacement of the beam in terms of the contact force, the forward Laplace transform is applied over time to the equation of motion of the computational domain, resulting in the following expression: where ˆ w c represents the unknown Laplace-domain displacement of the computational domain, and s = σ + i ω is the Laplace variable, with σ a small and positive real number and ω the angular frequency. Furthermore, ˆ where ˆ Q c represents the Laplace-domain contact force in the computational domain. It must be noted that the contact force ˆ Q c is yet unknown, while the initial-conditions forcing is completely determined ( Eq. (11) ). It is therefore chosen to determine the response caused by the moving load and the response generated by the initial conditions separately, using different approaches. The displacement ˆ w IC c generated by the initial conditions can be determined in the Laplace domain and the procedure is elaborated in Section 2.2.2 . The displacement caused by the unknown contact force can be obtained in the time domain by means of the Green's-function method. To this end, the Laplace-domain displacement ˆ w ML c caused by the unknown contact force is expressed in terms of the Green's function as follows: where g represents the time-domain Green's function and is obtained by applying the inverse Laplace transform to ˆ g , as shown in Section 2.2.2 .
The total displacement of the beam in the time domain is a superposition of the displacement caused by the moving load and the displacement generated by the initial conditions: where w L c represents the displacement of the computational domain caused by the oscillator continuing its movement in the right domain, and w L c ∼ H(v t − L ) meaning that w L c is non-zero only after the oscillator has entered the right domain. Because the oscillator is assumed to be in the steady state when it enters the right domain, w L c depends on the steady-state contact force Q 0 and not on the unknown contact force Q c (t) . Therefore, w L c is obtained using the same procedure employed for obtaining w IC c , which is elaborated in Section 2.2.2 . The displacement of the beam has now been expressed in terms of the unknown contact force, as seen in Eq. (19) . Taking advantage of the linearity of the vehicle model (except for the contact spring), the same procedure is applied for obtaining the displacement of the mass as a function of the unknown contact force. The Laplace-domain displacement of the mass ˆ u reads The time-domain displacement u of the mass is obtained by evaluating the inverse Laplace transform of Eq. (20) analytically, and its expression reads Both the displacement of the mass and that of the beam are now expressed in terms of the unknown contact force. The contact force Q c (t) is solved for iteratively at each time moment such that the contact relation is satisfied. The procedure for obtaining the contact force is thoroughly explained in [43] , and is briefly summarized here. Outside the computational domain, the contact force is assumed to be constant and equal to Q 0 . Therefore, the contact force needs to be determined only in the time interval when the oscillator is inside the computational domain. Upon evaluating the displacement of the beam at x = v t, both displacements ( Eqs. (19) and (21) ) are substituted in the contact relation ( Eq. (3) ), resulting in an integral equation: where H R (t) denotes the Heaviside function with argument R (t) . Considering the time discretization t 0 ,t 1 , ..., t n , with t 0 = 0 ,t n = t (the present time moment) and It must be emphasized that although Eqs. (22) and (23) are valid only for 0 ≤ t n ≤ L v , the simulation is not limited to this time interval provided that the oscillations of the mass have stopped before the oscillator reaches the right boundary x = L .
Assuming that both the Green's function of the beam-foundation system and the contact force have a linear variation during the time interval τ ∈ [ t i −1 , t i ] , the integrals in Eq. (23) are evaluated analytically, resulting in a nonlinear algebraic equation which can be solved iteratively for the unknown contact force at present time t n . Once the time history of the contact force has been obtained, the response of the mass and that of the beam can be computed by substituting the determined contact force in Eqs. (19) and (21) . At this point, the Green's function of the inhomogeneous and infinite beam-foundation sub-system g, the displacement w IC c caused by the initial conditions, and the displacement w L c generated by the oscillator continuing its movement in the right domain have not yet been determined. This is done in the next sub-section.

The beam-foundation Green's function, the response to initial conditions and the response caused by the oscillator continuing its movement in the right domain
Due to the spatial variation of the foundation stiffness and damping (as seen in Eq. (14) ), the Green's function ˆ g of the beam-foundation sub-system, the displacement ˆ w IC c caused by the initial conditions, and the displacement ˆ w L c generated by the oscillator continuing its movement in the right domain cannot be solved analytically for all stiffness and damping profiles. Therefore, the finite difference method (central difference scheme of O ( x 2 ) ) is used to spatially discretize Eq. (14) .
Upon discretization of the computational domain, the Laplace-domain equation of motion describing the beam-foundation sub-system subjected to initial conditions reads where ˆ K i j represents the bending-stiffness matrix of the beam, ˆ F D c ,i represents a force exerted by the semi-infinite domains at the boundaries derived in Section 2.2.3 , and I i j is the identity matrix. The time-domain displacement is obtained by left multiplying Eq. (24) with the inverse of the dynamic stiffness matrix (expression in the square brackets) and by evaluating the inverse Laplace transform numerically. However, the Laplace-domain spectrum of ˆ w IC c , j exhibits a poor decay due to the applied non-zero initial conditions. Consequently, the numerical integration must in principle be performed up to very high frequencies leading to a significant computational effort. A method of incorporating the high frequencies without increasing the computational effort is presented in [51] and is also applied here to accurately obtain w IC c , j . The same procedure used for obtaining w IC c is used to obtain w L c . The Laplace-domain equation of motion reads where ˆ F L c ,i represents a boundary-forcing vector exerted at the right boundary by the oscillator as it continues its movement in the right domain, and is derived in Section 2.2.3 .
To obtain the Laplace-domain Green's function ˆ g needed in Eq. (17) , the right-hand side of Eq. (14) is replaced by a Dirac delta function in space EI −1 δ(x − ξ ) . Upon spatial discretization, the Laplace-domain equation of motion that describes the impulse response reads The discretized Dirac delta function is a vector with a single non-zero entry at x = ξ , and its expression is given as follows: It must be noted that Eq. Firstly, in order to evaluate the Green's functions at x j = v t j and ξ i = v τ i ( Eq. (22) ) the sampling requirements x = v t and ξ = v τ need to be satisfied. Furthermore, depending on the required maximum frequency of integration ω max , the time step can become very small leading to an unrealistically small spatial step. To fulfill the above-mentioned requirements, a spatial interpolation scheme is used for ˆ g i j . Secondly, for the total response w c to instantaneously reach the steady state at t = 0 , the Green's functions must satisfy trivial initial conditions at t = τ, especially close to t = 0 . To this end, the causal character of the response is employed and the imaginary-part approach is selected since it automatically satisfies the trivial initial conditions even if the inverse Laplace integral is truncated. Consequently, the time-domain Green's functions are computed as At this point, the only unknowns are the boundary conditions of the computational domain. Upon application of the finite difference scheme, the boundary conditions are incorporated in the stiffness matrix of the beam ˆ K i j and in the boundaryforcing vectors ˆ F D c ,i and ˆ F L c ,i . The non-reflective boundary conditions are determined in the next sub-section.

Non-reflective boundary conditions
The procedure in this sub-section is based on the formulation described thoroughly in [51] , and is summarized in the following. As boundary conditions for the computational domain, the reaction forces of the semi-infinite domains at the interfaces with the computational domain are imposed. The goal is to express the interface reaction forces (bending moment and shear force) of the left and right domains as functions of the unknown displacement and slope of the computational domain at the corresponding interfaces. To this end, the forward Laplace transform is applied over time to the equations of motion of the two semi-infinite domains, ( Eq. (1) for x ∈ (−∞ , 0] and x ∈ [ L, ∞ ) ), leading to where ˆ w l and ˆ w r represent the unknown Laplace-domain displacements of the left and right domains, respectively; k l and k r represent the wavenumbers of the two semi-infinite domains and read  15) and (16) , but with Q 0 instead of Q c x v and with the initial conditions of the left domain ( Eq. (10) ) instead of those of the computational domain. As for the boundary conditions, the zero-displacements condition (the Laplace-domain counterpart of Eq. (8) ) is imposed at infinity, while the unknown Laplace-domain displacement and slope of the computational domain (the Laplace-domain counterparts of Eqs. (4) and (5) ) are prescribed at the interfaces.
The Laplace-domain displacements ˆ w l and ˆ w r can be obtained analytically by solving Eqs. (29) with the above-discussed boundary conditions. The interface reaction forces of the two semi-infinite domains are then expressed as functions of the prescribed displacements and slopes, and they read where the entries of the matrices represent the dynamic stiffness coefficients associated with boundary forces proportional to the unknown displacement and slope at the interfaces; subscript V stands for shear force, M for bending moment, υ for translation, and φ for rotation. The coefficients are given by Eqs. (46) and (47) of [51] . In addition, vector D l incorporates the influence of the initial conditions on the reaction forces, giving rise to boundary forces independent of the unknown displacement and slope of the computational domain, which are accounted for in Eq. (24) (24) ) and the response ˆ w ML c caused by the moving oscillator are treated separately, vector D l is only incorporated in the boundary condition for the determination of the response caused by the initial conditions. The vector is given by the following expression: where ˆ w l , p (0 , s ) is the particular solution that accounts for ˆ F IC l in Eq. (29) . It can be obtained as follows: where ˆ g l is the Laplace-domain Green's function of a homogeneous and infinite beam-foundation system with the properties of the left domain, and its expression is given by Eq. (54) of [51] . The integral in Eq. (33) can be computed analytically, but the solution is not presented here for brevity. Furthermore, ˆ V L and ˆ M L in Eq. (31) are the shear force and bending moment, respectively, exerted by the moving oscillator on the right boundary after it has entered the right semi-infinite domain, and their expressions read It must be noted that Eq. (34) is correct only if the oscillations of the vehicle have stopped before reaching x = L (as assumed in Section 2.1 and explained in Section 2.2.1 ), such that the contact force is again constant and equal to Q 0 ; this represents the criterion for choosing L . Also, Eq. (34) is used in the boundary conditions only for determining the displacement ˆ w L c generated by the oscillator continuing its movement in the right domain, and not for determining ˆ w IC c , nor for determining ˆ g . It must also be mentioned that for the results presented in Section 3 , Eq. (34) is not used because the simulations are performed only until the moving oscillator reaches the right boundary. However, according to the problem statement in Section 2.1 , if a certain analysis is concerned with the response of the computational domain after the moving oscillator has entered the right domain (i.e., free vibration), Eq. (34) must be incorporated. Moreover, the presented solution method can handle unstable vibrations of the oscillator onset by the inhomogeneity (e.g., breaking the wave-velocity barrier when going from the stiff domain to the soft one). However, in this case the maximum time of the simulation is limited to the moment the load reaches the right boundary of the computational domain because Eq. (34) is only valid for a constant contact force.
Eqs. (30) and (31) represent the non-reflective boundary conditions which are prescribed to the computational domain through the Laplace-domain image of Eqs. (6) and (7) . The term D l is accounted for through the boundary-forcing vector ˆ F D c ,i whileˆ V L andˆ M L are accounted for through the boundary-forcing vectorˆ F L c ,i , because they are not proportional to the unknown displacement and are thus considered as external forces. The remaining parts of Eqs. (30) and (31) are accounted for in the beam's stiffness matrix ˆ K i j . The application of these boundary conditions ensures that the behaviour of the finite computational domain is not influenced by artificial boundaries and that the infinite extent of the model is respected in an exact manner.

Results and discussion
In this section, the proposed solution is firstly validated by considering two limit cases and comparing the obtained results to a benchmark solution. Then, the time-domain response of the wheel and of the rail under the wheel are presented for two specific cases, and the influence of the oscillator on the response at the contact point is highlighted. Subsequently, the influence of the oscillator on the power input into the beam-foundation system is discussed. From these results, two indicators of potential damage caused to the supporting structure are identified, namely the maximum contact force and the energy/power input. Finally, the influence of the transition length, oscillator velocity, and stiffness ratio on the two indicators is assessed through a parametric study.
To investigate the influence of the transition length l t on the maximum contact force and the energy/power input, the spatial profile of the foundation stiffness is chosen based on a sine squared (depicted in Fig. 2 ) as follows: where x tc represents the location of the centre of the transition and p is the stiffness ratio between the stiff and soft domains. The location of the transition zone inside the computational domain is dictated by the spatial extent of the initial conditions, as discussed in Section 2.1 , while the value of L is chosen such that the oscillations of the moving mass have stopped before it reaches the position x = L, as discussed in Section 2.2.3 . In addition, the spatial variation of the foundation damping is chosen such that the damping ratio ζ , defined similarly to that of a single-degree-of-freedom system, is kept constant. Consequently, the foundation damping profile reads

Parameter values
The choice of two parameters, namely the stiffness ratio p and the velocity v of the oscillator, requires special attention.
In this work, the stiffness ratio p = is the stiffness ratio between the stiff and soft domains of the supporting structure (i.e., not including the beam's stiffness). This is different from the dissimilarity q of the vertical stiffness of the track, which is the difference in stiffness experienced by the wheel (i.e., including the beam's stiffness). The two are related (in static or quasi-static conditions) through the following expression [52] : Several researchers have measured the vertical track stiffness in areas with soft soils and have reported differences in the vertical track stiffness at transition zones between q ≈ 1 . 5 -3 [12,14,53,54] , which corresponds to p ≈ 1 . 5 -4.5. For most of the results in the remainder of the paper, p = 5 is chosen. Although it is on the high end, it is not unrealistically high. The value was chosen to accentuate the transition radiation mechanisms such that they are easier to identify. In Section 3.3.1 , p = 10 (i.e., q = 5 . 6 ) is chosen to illustrate the wheel-rail separation. Moreover, in Section 3.5.3 , p = 1 -20 (i.e., q = 1 -9.5); although from p = 5 onwards the values are less realistic (or even not at all realistic close to p = 20 ), they are presented in order to observe the overall trend. As for the velocity of the oscillator, its value is presented in this paper as a fraction of the critical velocity (i.e., minimum phase velocity c ph , min for this 1-D system). The Kelvin foundation is known to significantly overestimate the critical velocity of the railway track; for example, for the parameters used in this paper, c ph , min ≈ 410 m/s (i.e., 1475 km/h). This  is much higher than the measured critical velocity of railway tracks in areas with soft soils, where the so-called Rayleigh wave velocity (which usually is the critical velocity) can be as low as 40-50 m/s (i.e., 150-180 km/h) [55,56] . Therefore, choosing a realistic value for the oscillator velocity (e.g., 10 0-20 0 km/h) together with the Kelvin model that overestimates the critical velocity will lead to a quasi-static response equivalent to the train velocity of 20 km/h on a track with realistic minimum critical velocity of 180 km/h. Consequently, in this paper, the oscillator velocity is chosen relative to the c ph , min to represent realistic behaviour. A viable alternative is to incorporate the rail, sleepers, and ballast into the beam element with equivalent mass, bending stiffness, and internal damping [57] . Doing so, the critical velocity will have more realistic values; however, this model has its own shortcomings (e.g., obtaining the equivalent bending stiffness and internal damping is not straightforward). Nonetheless, this will not change the qualitative observations since, as explained above, the ratio of the train's velocity to the critical one is important, and not so much the absolute values of these quantities. The maximum design velocity of trains in areas with soft soils is around 60-70% of the minimum critical velocity in the track. However, in extreme cases, trains have even exceeded the critical velocity (e.g., the Swedish X-20 0 0 high-speed train which runs along the West Coast Line between Gothenburg and Malmö [56] ). In this paper, for Sections 3.3 and 3.4 , the oscillator velocity is chosen as 90-95% of c ph , min ; although these velocities are high for regular train operation, they are chosen to once more accentuate the transition radiation mechanisms such that they are easier to identify. It must be emphasized that the identified mechanisms are present also for lower velocities, but they are just harder to observe. For Section 3.5 , the oscillator velocity is varied between 13-110% of c ph , min . Again, although super-critical velocities are rarely observed in practice, the wide range of velocities is given to analyze the overall trend in behaviour.
The parameters which are kept constant throughout the presented results are given in Table 1 . The bending stiffness corresponds to the UIC 60 rail and the mass per unit length includes the rail mass per meter and the contribution of the (concrete) sleepers uniformly distributed along the rail [58] , such that the model captures the track response at low frequencies (including its first natural frequency).

Validation and convergence
To validate the solution derived in Section 2 , two limit cases are considered and are compared to a benchmark solution. Firstly, to validate the Green's functions of the beam-foundation structure and the suppression of initial transients through the initial conditions applied, the limit case of a moving constant load is considered, where the magnitude of the load is equal to Q 0 . In this limit case, the degree of freedom of the mass and the contact equation vanish, while the time-domain displacement of the beam ( Eq. (19) ) simplifies to Secondly, to incorporate the iterative scheme of solving the contact equation in the validation, the same limit case is considered (moving constant load) by taking the limit of the wheel's mass going to zero. However, the wheel's mass cannot be set to zero because it leads to a singularity in the expression of the wheel's displacement, as can be seen in Eq. (21) . Therefore, the wheel's mass is gradually decreased to study the convergence of the second limit-case solution w B lim as the mass goes to zero. The foundation is considered as piecewise homogeneous for both cases.
The limit-case responses are compared to the semi-analytical solution for the piecewise homogeneous system subjected to a moving constant load. This benchmark solution is obtained by solving a system of two semi-infinite domains with different foundation stiffness by using the Fourier Transform over time. In the Fourier domain, the displacements of the two domains can be obtained analytically by imposing interface conditions and the condition of zero displacements at infinity [e.g., 22,50 , 51 ]. To obtain the time-domain solution, the inverse Fourier integral is evaluated numerically.
The errors e A and e B presented in Fig. 3 , corresponding to the first and second limit cases, respectively, are defined as follows: where w bench is the benchmark solution. To study the convergence of the presented solution, the truncation frequency for the inverse Laplace transform in the first limit case is varied, while a constant truncation frequency of f max = 16 kHz is used for the benchmark solution. Note that by changing the maximum frequency for the limit case, according to the Nyquist sampling rule, the time stepping also changes. Panel (a) of Fig. 3 shows that the first limit-case solution converges to the semi-analytical one as the maximum frequency is increased. The error in the soft domain is smaller than the one in the stiff domain because the higher foundation stiffness leads to smaller displacements and, thus, to a higher relative error. The error e A is mainly caused by two factors, namely the truncation (in frequency) of the inverse Laplace integral ( Eq. (28) ) and the spatial discretization introduced by the finite difference method. Also, to satisfy x = v t (discussed in Section 2.2.2 ), the spatial step size x can become unrealistically small depending on the time step t. To this end, a spatial step size which leads to accurate results is chosen ( x = 0 . 05 m), and to satisfy x = v t an interpolation scheme is used. However, the error caused by the interpolation scheme is negligible. It can be observed in Fig. 3 that for truncation frequency f max = 8 kHz, the error is less than 0.5%, which is more than satisfactory. A higher maximum frequency leads to smaller error; however, the computational effort increases significantly. For all other results presented in this section, the maximum frequency was therefore chosen as 8 kHz.
In panel (b) of Fig. 3 , the second limit-case solution w B lim is observed to converge to the benchmark solution for the wheel's mass tending to zero. Moreover, the errors in the steady-state regimes are the same for all mass values and are equal to the errors in the first limit case e A . This suggests that the error introduced by iterative solver is negligible. It must be noted that error e B in the area of the transition zone is mainly due to the incapability of considering the true limit case ( M = 0 ), and does not imply anything about the error in the cases studied further, where the inertia of the wheel is significant.
To investigate if the whole solution is valid and not just at the contact point, the full displacement field in the limit case A (i.e., Eq. (38) , but not limited to x = v t) is presented for different time moments in Fig. 4 , and is compared to the benchmark solution w bench (x, t ) . It can be seen from Fig. 4 that the results agree very well from the fact that the two lines are indistinguishable. In panel (a) it can be seen that the initial transients have been successfully suppressed through the implementation of the correct initial conditions (i.e., the response is in the steady-state regime); here it can also be clearly observed why the initial conditions based on the eigenfield must decay before reaching the inhomogeneity ( x tc = 10 m in this case). Moreover, in panel (e), the fact that the two solutions agree very well demonstrates that all waves generated inside the computational domain propagate away from the transition with no reflection at the boundaries. Also, in panel (f) it can be seen that the eigenfield passes to the right domain without any disturbance. Therefore, it can be concluded that the solution derived in Section 2 is correct. Finally, a validation of the vehicle model employed in this paper is presented in Fig. 5 . More specifically, the suitability of the loaded one-mass oscillator for studying the wheel-rail interaction is investigated. The normalized maximum contact force obtained with the loaded one-mass oscillator is compared to the one obtained with a more comprehensive vehicle model, namely the three-mass oscillator used by Ang and Dai [19] . For the higher velocity (panel (a)), the results agree very well, while for the lower velocity (panel (b)), the results present some differences, but not major ones. This can be explained by the fact that the lower the velocity, the lower the frequencies of interaction, which leads to the sprung masses having a more significant influence. Overall, the results agree well rendering the loaded one-mass oscillator suitable for studying wheel-rail interaction. Nonetheless, it must be noted that a vehicle model with multiple contact points can influence the results more significantly. Also, for large transition lengths, the inertia of the suspended masses may have a more pronounced influence due to the low frequencies of interaction.

Time-domain responses
To study the interaction between the moving oscillator and the beam-foundation structure, the displacement of the mass and that of the beam at the contact point, as well as the contact force are presented in Fig. 6 . Two specific cases are considered: the oscillator moving from the soft domain to the stiff one and that moving from the stiff domain to the soft one. For both cases, soft-to-stiff and stiff-to-soft, the velocity is chosen as 90% of the minimum phase velocity in the soft domain ( c l , ph , min for the soft-to-stiff case and c r , ph , min for the stiff-to-soft case), and the transition length l t is chosen as   [19] with the three-mass oscillator (circles) for v = 90 m/s (panel (a)) and v = 70 m/s (panel(b)); the parameters used are the ones from Ang and Dai [19] . 0.1 m, which is close to the piecewise-homogeneous case. To ensure that the initial displacement and velocity fields do not interact with the inhomogeneity, the centre of the transition zone x tc is positioned at x = 19 m for the soft-to-stiff case and at x = 11 m for the stiff-to-soft case. It is clear from Fig. 6 that the response reaches the steady state instantaneously as the oscillator enters the computational domain and thus, that the initial transients have been suppressed. Note that this initial plateau is only depicted in Fig. 6 to demonstrate that the initial transients have been suppressed, and is omitted in the following figures since it does not provide any information about the interaction.
In the soft-to-stiff case (left panels in Fig. 6 ), the displacement of both the mass and the beam (at the contact point) exhibit significant oscillations (panel (a) in Fig. 6 ). The oscillations are caused by the interaction of the moving oscillator with the inhomogeneous elastic structure, where three phenomena play a role: 1. Bending waves are excited as the load approaches and passes the transition due to the eigenfield interacting with the inhomogeneity. 2. A vibration of the mass (i.e., variation of the vertical momentum) induced by bending waves that kinematically excite the moving oscillator at the contact point. 3. A vibration of the mass induced by a parametric variation of the system properties as experienced by the moving oscillator (i.e., varying foundation stiffness).
The three phenomena described here are coupled as follows. The first and second phenomena are completely interdependent; stronger wave radiation introduces more oscillations at the contact point, leading to a stronger vibration of the mass, which in turn influences the wave radiation. In addition, although stronger vibration of the mass leads to stronger wave radiation, implying that the third phenomenon influences the first one, stronger wave radiation does not lead to stronger parametric variation, meaning that the first phenomenon does not influence the third one. Furthermore, the first phenomenon is considered to be the primary one, as it also takes place in the case of a moving constant load [51] , while the second and third phenomena are characteristic of the moving oscillator case.
The stronger radiation caused by the oscillation of the mass has a clear influence on the variation of the contact force (panel (c) in Fig. 6 ). In the considered case, the maximum contact force is almost three times larger than in the steady state. Moreover, the maximum peak in the contact force is followed by a decrease in contact force which in this case reaches the value of zero for a short time interval. This means that the wheel loses contact with the rail. The loss of contact between the wheel and the rail is discussed more thoroughly in Section 3.3.1 . The large variation in the contact force could be one of the factors leading to degradation of the foundation in the vicinity of transition zones. Therefore, the maximum contact force is selected in Section 3.5 as an indicator of the damage in the supporting structure and thus of the performance of the transition zone.
In the stiff-to-soft case (right panels in Fig. 6 ), the displacement of both the mass and the beam (at the contact point) exhibit smaller oscillations compared to the soft-to-stiff case, which is in agreement with the literature [19,51] . This can be explained based on the phenomena described previously: a) The smaller in displacement and less broad eigenfield turns out to yield weaker bending-wave excitation at the inhomogeneity and the excited large-amplitude waves propagate into the soft domain (related to the first phenomenon). b) Smaller oscillations at the contact point lead to smaller vibrations of the mass (related to the second phenomenon), which in turn introduce weaker additional waves in the beam (related to the first phenomenon). c) As the oscillator passes the transition, the displacement at the contact point increases due to lower foundation stiffness; however, the inertia of the mass keeps it (initially) at the same level. Consequently, in the stiff-to-soft case, the inertia force points upwards while the dead load points downwards, partially cancelling each other (related to the third phenomenon).
The effect of items b) and c) are clearly visible in the lower variation of the contact force in the stiff-to-soft case (panel (d) in Fig. 6 ).
It is important to note that the stiff-to-soft case analyzed in this paper assumes that the system initially is in the steady state. In reality, the system reaches the steady state in the stiff domain (e.g., a bridge) only if the stiff domain is long enough. For a short stiff domain embedded in a soft domain (e.g., a culvert [30] ), the response might not reach the steady state, leading to different behaviour in the stiff-to-soft transition.
To highlight the influence of the moving oscillator on the vibrations at the contact point, the displacement of the beam at the contact point is compared to the displacement of the beam under a moving constant load in Fig. 7 . The steady-state responses are the same in both cases (oscillator vs. constant load) because the vertical acceleration of mass is zero in this regime. The differences arise once the eigenfields start to interact with the inhomogeneity. In the soft-to-stiff case (panel (a) of Fig. 7 ), due to the decrease in contact force (panel (c) in Fig. 6 ), a smaller absolute displacement is observed at the contact point before the transition ( x ≈ 18 m). The contact force is lower because the reflected wave destructively interferes with the eigenfield, thus unloading the contact between the wheel and the rail. Here, the oscillator can have a beneficial influence on the extreme displacements. However, once the moving vehicle has entered the right domain, significant vibrations are observed in the oscillator case. Here, the oscillator has a negative influence on the extreme displacements. Moreover, the larger oscillations could lead to more significant rearrangements of the ballast particles, potentially leading to differential settlements. In the stiff-to-soft case, the difference between the moving oscillator and the moving constant load is significantly smaller. This could have been anticipated from the smaller variation in the contact force presented in panel (d) of Fig. 6 .

Loss of contact between the wheel and the rail
For certain values of the oscillator velocity, stiffness ratio and transition length, the wheel can loose contact with the rail. This phenomenon has been identified in studies of transition radiation [e.g., 19] , of oscillators interacting with bridges [e.g., 59] , and in studies of stability of moving vehicles [e.g., 45] . The loss of contact occurs in the case considered in the previous section (soft-to-stiff transition), but for a limited duration. To highlight it, the stiffness ratio is increased to p = 10 for the results presented in this sub-section (3.3.1) .
The left panels in Fig. 8 clearly present the loss of contact between the wheel and the rail, which occurs when the contact force (panel (c) in Fig. 8 ) is zero. In the considered case, there are even two loss-of-contact events. This can also be observed in the displacements (panel (a) in Fig. 8 ) and velocities (panel (e) in Fig. 8 ). When the displacement and velocity of the beam at the contact point have similar patterns as the displacement and velocity of the mass, respectively, there is contact between the two. Once the contact is lost, the behaviour of the beam and the mass become considerably different. This is due to the different forces acting on the two sub-systems when there is loss of contact; the beam is subject to the restoring forces of the foundation and of the bending resistance of the beam, while the mass is acted upon by the dead load Q 0 . The loss-of-contact phenomenon and the impact after it lead to strong oscillations in the response and large variation in the contact force, potentially causing damage in the supporting structure and/or in the rail.
The right panels in Fig. 8 show that the loss of contact between the wheel and the rail does not occur in the stiff-to-soft case (for the chosen parameters). This can be explained by the same reasoning as employed in the previous sub-section.

Energy considerations
In this section, the vehicle-structure interaction is studied from the energy point of view, which offers additional insight into the influence of the oscillator on the transition radiation and thus into the potential damage in the supporting structure. The energy emitted by the vehicle into the track that is dissipated in the ballast layer can be directly related to the damage in the supporting structure [60] . The higher the energy emitted (i.e., energy input) into the track, the stronger the degradation. Consequently, reducing the energy input represents an important aim for damage reduction in the supporting structure. The current model cannot be used to identify the amount of energy dissipated in the ballast layer. However, the proposed model can provide the energy input into the track, which can be an estimate of the energy dissipated in the ballast layer, and thus an indicator of potential damage in the supporting structure. To highlight the influence of the oscillator on the energy input, the oscillator is compared to the moving constant load.
The energy input per unit time, also referred to as the power input, provided by the moving oscillator P i , osc and by the moving constant load P i , const are given by the following expressions: where w 0 is the displacement of the beam in the moving constant load case. The total energy spent by the moving oscillator and the moving constant load can be obtained by integrating the power inputs in Eq. (40) over the whole time domain. However, these integrals are not convergent due to the damping in the foundation, which requires continuous power input to maintain the constant velocity. For that reason, the differential power input between the moving oscillator and the moving constant load can be integrated instead (done below and in Section 3.5 ), to offer information about the difference in energy input. In addition, one can look at the spectrum of the differential-power input, which gives insight into the frequency content of the additional energy exchange between the oscillator and the supporting structure. The spectrum ˜ P i can be computed by taking the forward Fourier transform of the differential-power input, as given in the following expression: Fig. 9 presents the power inputs for the two vehicle models and the amplitude spectra of the differential-power input. Since the inertia of the mass is not activated in the steady-state regime, the power inputs for the two vehicle models are initially the same. When the eigenfields interact with the inhomogeneity, the power inputs fluctuate significantly before reaching the steady-state regime in the right domain. The steady-state power input for the soft foundation is higher than that for the stiff foundation because the slope of the rail under the moving oscillator (moving load) is larger as the foundation is softer. Overall, the power input in the soft-to-stiff case appears to be significantly amplified by the presence of the oscillator, thus showing that the presence of the oscillator can lead to stronger transition radiation. For the stiff-to-soft scenario, the influence of the oscillator is less drastic.
As can be seen from the amplitude spectra of the differential-power input (bottom panels of Fig. 9 ), a significant part of the additional energy exchange in the soft-to-stiff case is at frequencies larger than both cut-off frequencies (of the soft and stiff domains). For the stiff-to-soft scenario, the additional energy exchange lies in similar amounts below and above the cut-off frequency of the soft domain. To give a reference, the amplitude spectrum of the eigenfield's power flux in the left domain is included in the bottom panels of Fig. 9 . It can be seen that for the soft-to-stiff case, the additional power exchange between the oscillator and the supporting structure is significantly larger than the spectrum of the eigenfield's power flux while for the stiff-to-soft case they are of similar magnitude.

Parametric study
In Sections 3.3 and 3.4 , two indicators of potential damage caused to the supporting structure have been identified, namely the maximum contact force and the difference in energy input. In this section, these two indicators are addressed as functions of the oscillator velocity, the transition length and the stiffness ratio through a parametric study. The three  parameters chosen to be varied are the most influential ones for the transition radiation phenomenon. Furthermore, these parameters can be adjusted in the design stage to minimize damage in the supporting structure. It must be noted that, while the maximum contact force offers insight into the magnitude of the transition radiation (as a whole), the difference in energy input provides information only into the additional effect of the vehicle-structure interaction on the transition radiation.

Variation of the transition length
Here, the transition length l t is varied and the response for different oscillator velocities is investigated in terms of the maximum contact force ( Q max = max (Q (t)) ) and in terms of the difference in energy input. It must be emphasized that in the present study the transition length refers to the zone over which the stiffness (and damping) of the supporting structure varies, and that the rail is considered completely straight (i.e., no differential settlements) when it is at rest, while in some other studies (e.g., Lei and Mao [26] ) the transition length/angle refers to the zone where the rail exhibits some slope due to differential settlements. Fig. 10 presents the maximum contact force versus transition length for different oscillator velocities. In the soft-to-stiff case (panel (a) in Fig. 10 ) the maximum contact force decreases nonlinearly with increasing transition length for all velocities considered. The decrease is most significant for the velocities close to the minimum phase velocity. However, even for operational speeds of conventional trains (e.g., v = 0 . 5 c ph , min ), the maximum contact force can be significantly minimized Fig. 11. Difference in energy input versus transition length for different oscillator velocities; the soft-to-stiff case (panel (a)) and the stiff-to-soft case (panel (b)); in the legend, c ph , min represents c l , ph , min for the soft-to-stiff case and c r , ph , min for the stiff-to-soft case; p = 5 .
by choosing a transition length of 4-6 m. Increasing the transition length even more does not lead to a considerable improvement of the transition zone performance. In the stiff-to-soft case (panel (b) in Fig. 10 ), the maximum contact force decreases as well with increasing transition length, but much less pronounced than in the stiff-to-soft case. Nonetheless, for small transition lengths, the maximum contact force can be 20-25% larger than the steady-state one for relatively high velocities, and the maximum contact force can be reduced by 10-15% when increasing the transition length to 6-10 m. Fig. 11 presents the difference in energy input between the moving oscillator and the moving constant load versus transition length for different velocities. A similar behaviour as for the maximum contact force is observed in the soft-to-stiff case (panel (a) in Fig. 11 ). The difference in the energy input decreases nonlinearly with increasing transition length and the decrease is drastic for relatively high velocities, about a factor 2-5 in going from l t = 0 . 1 m to l t = 10 m. For velocities up to v = 0 . 75 c ph , min , the behaviour in the stiff-to-soft case (panel (b) in Fig. 11 ) is similar to the soft-to-stiff case. However, for larger velocities (e.g., v = 0 . 95 c ph , min and v = 1 . 10 c ph , min ), the difference in the energy input exhibits a behaviour totally different from the soft-to-stiff case and from the maximum contact force. A peak in the difference of energy input is identified between l t = 2 m and l t = 6 m. This means that in this range, although most can be gained in terms of the maximum contact force decrease (see panel (b) of Fig. 10 ), the energy input increases compared to the constant load case. This could be one of the reasons why similar amount of damage is observed in front and behind a transition zone consisting of a softto-stiff transition directly followed by a stiff-to-soft one (for one-way tracks), where mitigation measures to create a gradual transition are employed (e.g., approach slabs). In the case of an approach slab of 2-4 m, there is a substantial decrease in the difference of energy input before the transition (compared to the abrupt case), while there is a substantial increase in the difference of energy input after the transition. This finding is not obvious from studying solely the observables such as contact force, displacements and velocities.

Variation of the velocity
In this section, the maximum contact force and the difference in energy input are presented as functions of the velocity of the moving oscillator. Three transition lengths are considered: l t = 0 . 1 m, l t = 5 m and l t = 10 m. Fig. 12 presents the maximum contact force results. In the soft-to-stiff case (panel (a) in Fig. 12 ), the maximum contact force increases nonlinearly with the increasing velocity up to a value close to the minimum phase velocity c l , ph , min , and it decreases beyond that. For an abrupt transition, the increase of the maximum contact force is significant, but it reduces with increasing transition length. This is in accordance with the findings in Section 3.5.1 ( Fig. 10 ). Furthermore, the critical velocity for the maximum contact force (marked with rectangles in Fig. 12 ) is observed to decrease with increasing transition length. For l t = 0 . 1 m the critical velocity is close to c l , ph , min , while in the case of l t = 10 m the critical velocity has decreased to 0.9 c l , ph , min . A different behaviour is observed for the stiff-to-soft cases (panel (b) in Fig. 12 ). A local peak in the maximum contact force is observed at relatively low velocities, after which an almost linear increase is observed with increasing load velocity. The location of the local peak shifts to higher velocities for increasing transition length. Also, compared to the soft-to-stiff case, the change in maximum contact force with increasing load velocity is much smaller.
The difference in energy input for the soft-to-stiff case (panel (a) of Fig. 13 ) exhibits a similar behaviour to the one observed for the maximum contact force (panel (a) in Fig. 12 ). However, the critical velocity for the difference of energy input shifts much less with the variation of the transition length. Moreover, the difference in energy input becomes negative at certain supercritical velocities, meaning that the moving oscillator inputs less energy into the beam-foundation system than the moving constant load. In the stiff-to-soft case (panel (b) of Fig. 13 ), the difference in energy input exhibits a behaviour different from the one observed for the maximum contact force (panel (b) in Fig. 12 ). The difference in energy input increases with increasing transition length, with a peak around c r , ph , min , beyond which it decreases. The critical velocity exhibits some peculiar shifts with the variation of transition length. Moreover, the transition length increase does not seem to have such a considerable impact on the reduction of the difference in energy input compared to the soft-to-stiff case.    (panel (a)) and for the stiff-to-soft case (panel (b)); l t = 0 . 1 m; in the legend, c ph , min represents c l , ph , min for the soft-to-stiff case and c r , ph , min for the stiff-to-soft case.

Variation of the stiffness ratio
In this section, the maximum contact force and the difference in energy input are presented as functions of the stiffness ratio p. Lei and Mao [26] found that the influence of the stiffness ratio on the maximum contact force is minimal. In their study, the supporting structure was modelled as two Kelvin-foundation layers in series; the bottom layer represents the subsoil and is piecewise homogeneous while the top layer represents the ballast and is homogeneous. Although the stiffness jump prescribed in the bottom layer is significant ( n = 10 and n = 100 , where n = k stiff lower /k soft lower ), the equivalent single-layer Kelvin foundation has a small stiffness ratio and a limited variation range ( p ≈ 1 . 6 for n = 10 and p ≈ 1 . 8 for n = 100 , while p = 5 can be observed in areas with soft soils [14,53,54] ). Moreover, Ang and Dai [19] found that for a smooth surface of the rail the increase in the maximum contact force with increasing stiffness ratio is negligible. However, the load velocities considered in their study are much lower than the minimum phase velocity in the soft domain, namely around 13% and 18%. Here, the velocity and the stiffness-ratio ranges are extended such that more general observations can be made about the influence of the stiffness ratio on the maximum contact force. Fig. 14 presents the maximum contact force versus stiffness ratio p for different oscillator velocities. In both cases (softto-stiff and stiff-to-soft), the maximum contact force increases with increasing stiffness ratio, tending to a specific value which is obtained in the limit case of k d , stiff going to infinity. Although the qualitative behaviour in the two cases is similar, Fig. 15. Difference in ener gy input between the moving oscillator and the moving constant load versus stiffness ratio for the soft-to-stiff case (panel (a)) and for the stiff-to-soft case (panel (b)); in the legend, c ph , min represents c l , ph , min (soft-to-stiff case) and c r , ph , min (stiff-to-soft case). the increase of the maximum contact force is much more significant in the soft-to-stiff case. In the soft-to-stiff case, the lines of v = 0 . 95 c ph , min and of v = 1 . 10 c ph , min intersect, implying that the critical velocity for the maximum contact force is also dependent on the stiffness ratio (next to transition length, discussed in Section 3.5.2 ). In the stiff-to-soft case, the line of v = 0 . 50 c ph , min lies above the lines of v = 0 . 75 c ph , min and of v = 0 . 95 c ph , min , which is in agreement with the results presented in panel (b) of Fig. 12 . Furthermore, Fig. 15 presents the difference in energy input versus stiffness ratio p for different velocities. The behaviour exhibited by the difference in energy input is very similar to the one observed for the maximum contact force ( Fig. 14 ).
As can be seen in the panel (a) of Fig. 14 , for small oscillator velocities (e.g, v = 0 . 13 c ph , min ), the increase in maximum contact force with increasing stiffness ratio is negligible, a result which is in agreement with the findings of Ang and Dai. (It must be mentioned that there is a difference between this agreement in findings and the one presented in Section 3.2 ; in Fig. 5 , the results are obtained using the same parameters as in Ang and Dai to validate the choice of the vehicle model, while here, the findings agree also for the parameter values adopted in the current paper.) Also, for stiffness ratio between 1.6-1.8, the change in maximum contact force is insignificant (less than 1% for small oscillator velocities and up to 4% for oscillator velocities close to the minimum phase velocity), results which are in agreement with the conclusions of Lei and Mao. However, for larger velocities and wider stiffness-ratio range, the maximum contact force increases drastically, reaching a dynamic amplification factor greater than 5 for velocities close theminimum phase velocity. Moreover, even for oscillator speeds of 25-50% of the critical velocity, which represent operational speeds of conventional trains, the dynamic amplification factor can reach values of 1.4-2.5 (for very large stiffness dissimilarities). Therefore, even if the rail is completely straight (i.e., no differential settlements have occurred), which is assumed in this paper, the stiffness variation alone can lead to a significant increase in the maximum contact force for large enough oscillator velocities and considerable stiffness-ratio values. This increase in the maximum contact force can lead to the onset of degradation in the form of differential settlements, which generally cause additional increase in the maximum contact force as has been found by Ang and Dai and Lei and Mao, aggravating the degradation. To conclude, these findings confirm that uneven stiffness, next to initial imperfection (as studied by others [e.g., 19 , 26] ), is one of the main causes for the onset of degradation.

Conclusions
In this paper, the influence of accounting for the interaction between the vehicle and the supporting structure on the transition radiation phenomenon has been studied. To this end, a one-dimensional model has been formulated, consisting of an infinite Euler-Bernoulli beam resting on a locally inhomogeneous Kelvin foundation, interacting with a moving loaded oscillator that has a nonlinear Hertzian spring. The linear behaviour of the individual sub-systems, namely the moving mass and the supporting structure, allows to obtain the solution by means of the Green's-function method. The finite difference method has been used for the spatial discretization of the finite computational domain to accommodate the smoothly inhomogeneous foundation. Furthermore, the infinite extent of the system has been accounted for through a set of nonreflective boundary conditions derived by replacing the semi-infinite domains adjacent to the computational domain with their response at the interfaces, and through appropriate initial conditions. The solution has been validated by comparing the obtained results to a benchmark solution for two limit cases in which the problem simplifies to a moving constant load, and by comparing the obtained results with findings in the literature.
Accounting for the interaction between the moving oscillator and the supporting structure generally leads to stronger transition radiation, caused by the variation of the vertical momentum of the moving mass. Results show that the interaction causes a variation of the contact force and the maximum contact force exhibits a significant increase compared to the moving constant load case. However, the maximum contact force decreases nonlinearly with increasing transition length. Furthermore, for oscillator velocities close to the minimum phase velocity, the maximum contact force increases significantly with increasing stiffness ratio, implying that the stiffness ratio is an important factor in minimizing the maximum contact force. This finding supplements the existing literature where limited velocity and stiffness-ratio ranges have been considered, while the present work provides a more complete picture about the influence of the velocity and stiffness ratio on the maximum contact force.
The influence of the oscillator on the transition radiation has also been investigated from an energy point of view. Results show that the power input by the oscillator into the supporting structure can be significantly larger than that of the moving constant load. Meanwhile, the difference in energy input between the moving oscillator and the moving constant load shows a significant decrease with increasing transition length for the soft-to-stiff transition. For the stiff-to-soft transition, the maximum difference in energy input does not necessarily occur for the minimum transition length, which is a peculiarity of this problem. Furthermore, the difference in energy input increases with increasing stiffness ratio following a similar trend to the one exhibited by the contact force.
Finally, two quantities, namely the maximum contact force and the difference in energy input, have been identified as possible indicators of damage occurring in the supporting structure. They are valuable quantities which can be computed in the preliminary design stages to assess the performance of a railway track transition zone.

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. It must be noted that the expressions presented in this appendix are not valid when Eq. (A.3) leads to repeated roots. This situation occurs when there is no viscous damping and v = c ph , min , or for very high viscous damping. These scenarios are not of relevance in this paper.