Transition radiation in a piecewise-linear and infinite one-dimensional structure–a Laplace transform method

Transition zones in railway tracks are areas with considerable variation of track properties (i.e., foundation stiffness), encountered near structures such as bridges and tunnels. Due to strong amplification of the response, transition zones require frequent maintenance. To better understand the underlying degradation mechanisms, a one-dimensional model is formulated, consisting of an infinite Euler–Bernoulli beam resting on a locally inhomogeneous and nonlinear Winkler foundation, subjected to a moving load. The nonlinearity is characterized by a piecewise-linear stiffness, and the system thus behaves piecewise linearly. Therefore, the solution can be obtained by sequentially applying the Laplace transform combined with the Finite Difference Method for the spatial discretization and derived non-reflective boundary conditions. Results show that the plastic deformation is a consequence of constructive interference of the excited waves and the response to the load’s deadweight, particularly for the soft-to-stiff transition. The plastic deformation area decreases quasi-monotonically with increasing transition length, and for super-critical velocities, small transition lengths and/or large stiffness dissimilarities, parts of the foundation experience plastic deformation even in the second loading–unloading cycle. Furthermore, the nonlinearity causes the maximum energy associated with the waves radiated forward and the maximum energy input not to occur for the smallest transition length, contrary to findings in corresponding linear systems. Moreover, the energy input drastically increases for the second passage of the moving load, making it a possible indicator of the damage in the supporting structure. The novelty of the current work lies in the computationally efficient solution method for an infinite system which locally exhibits nonlinear behaviour and in the study into the influence of the foundation’s nonlinear behaviour on the generated waves (i.e., transition radiation), and on the resulting plastic deformation. The model presented here can be used for the preliminary design of transition zones in railway tracks.


Introduction
Transition radiation is emitted when a source moves along a straight line with constant velocity and acts on or near an inhomogeneous medium [1,2]. It occurs, for example, when a train crosses an inhomogeneity in a railway track, such as a variation in foundation stiffness. The velocity of high-speed trains or conventional trains running on soft soils may approach the critical velocity in the supporting structure leading to strong transition radiation. Apart from potentially giving rise to vehicle instability and passenger discomfort, transition radiation has been addressed as one of the causes of track and foundation degradation due to the often strong amplification of the stress and strain fields [3][4][5][6]. This leads to a high frequency of maintenance required for transition zones in areas with soft soils, which can be 3-8 times higher than for the regular parts of the railway track [7,8]. Frequent maintenance leads to high costs and reduced availability of the track.
The first study on transition radiation of elastic waves was published by Vesnitskii and Metrikin [9]. It considered an infinite string on a piecewisehomogeneous Winkler foundation subjected to a moving point force and a moving point mass, respectively. To also account for the flexural rigidity of the system, Vesnitskii and Metrikin [10] analysed transition radiation in a semi-infinite beam on elastic springs clamped at one end. This system is equivalent to an infinite one where part of the foundation has infinitely high stiffness. Later on, the problem of a beam resting on a piecewise-homogeneous Winkler foundation was also addressed by others using different solution methods: modal expansion techniques [11,12] and the moving element method [13]. To study wave propagation in the ground due to transition radiation, 2-D models of a continuum acted upon by a moving load were analysed. The continuum was modelled as a piecewisehomogeneous half-space [14,15] and as a piecewisehomogeneous finite-depth layer [16,17].
The geometry of the transition zone has a strong influence on the transition radiation, and modelling it as piecewise homogeneous usually does not suffice for design purposes. The detailed geometry of the transition zone has been more accurately modelled using 2-D and 3-D finite-element models or combined boundaryelement models and finite-element models [e.g., [18][19][20][21][22][23][24]; a good review of models for transition zones can be found in [7]. Such models are necessary at the final stage of designing a transition zone, but they are computationally demanding. A simplified way to model a transition zone in a 1-D model, but more accurately than assuming a piecewise-homogeneous foundation, is to incorporate a smooth transition between domains with constant support stiffness. This has been introduced by Metrikine et al. [25] for the model of an infinite stringfoundation system. Incorporating a smooth transition in a 1-D model opens the possibility of studying the performance of different transition-zone geometries in earlier design stages.
Many studies about the dynamic response of railway tracks assume linear behaviour of the supporting structure. However, the influence of the foundations nonlinearity on a railway track is significant, as shown experimentally by Dahlberg [26], and should not be overlooked. A study of the steady-state response of an infinite string supported by nonlinear elastic springs subjected to two moving point loads can be found in [27]. To also account for the flexural rigidity, the steadystate response of a beam on a homogeneous and nonlinear elastic foundation subjected to a moving load has been analysed, considering a finite system [28][29][30] and an infinite one [31][32][33]. In addition, Hoang et al. [34] studied the steady-state response of an infinite beam with periodic nonlinear elastic supports.
Transition radiation has also been addressed in systems with nonlinear elastic behaviour of the supporting structure. Castro Jorge et al. [35] used a nonlinear elastic foundation to analyse the effect of the nonlinearity on the maximum displacements in a finite and piecewise-homogeneous system. In addition, Varandas et al. [36,37] considered nonlinear elastic behaviour of the supporting structure in a 3-D finite-element model of a transition zone. However, to study the degradation in a transition zone, the employed model should incorporate plastic behaviour of the foundation. To this end, Varandas et al. [38] developed a finite 1-D model describing the accumulated permanent deformation in a transition zone using a phenomenological model for the cyclic degradation of the supports. Moreover, Gallego Giner et al. [39] used an elastic-perfectly plastic model (i.e., Drucker-Prager) for the supporting structure in his study of a transition zone using a 3-D finite element model. Detailed 3-D models of finite and inhomogeneous systems that incorporate nonlinear behaviour of the foundation are available in the literature, as shown previously. However, simplified models of transition zones in infinite systems with nonlinear elasto-plastic foundation behaviour are not available in the literature. This motivates the aim of the current work, which is to formulate a 1-D model of an infinite Euler-Bernoulli beam on a smoothly inhomogeneous and nonlinear elastoplastic Winkler foundation, subjected to a moving load, and to study the effect of the nonlinear behaviour on the transition radiation and the degradation in the transition zone. The novelty of the current work is twofold.
Firstly, a computationally efficient solution method for an infinite system which locally exhibits nonlinear behaviour is presented. Secondly, the influence of the foundation's nonlinear behaviour on the generated waves (i.e., transition radiation), and the resulting plastic deformation are studied. The model presented here can be used for preliminary designs of transition zones in railway tracks by assessing the potential damage (i.e., plastic deformation) occurring in the transition zone as a function of the smoothness of the transition (i.e., length of the transition), of the moving load velocity, of the system's damping and of the stiffness dissimilarity. Furthermore, the results presented here offer insight into the physical mechanisms leading to degradation in the supporting structure.
This paper is structured as follows. In Sect. 2, the problem statement is presented and the solution is derived. Section 3 introduces the constitutive relation of the foundation stiffness and the assumptions considered. In Sect. 4, the results are presented and discussed, while Sect. 5 presents the drawn conclusions.

Problem statement
In this section, a 1-D model is formulated, consisting of an infinite Euler-Bernoulli beam resting on a smoothly inhomogeneous and nonlinear Winkler foundation, subjected to a constant moving load (Fig. 1). In view of practical relevance (degradation is often encountered inside or close to the transition zone), the nonlinear behaviour of the supporting structure is restricted to the transition zone and its vicinity; this domain is referred to as the computational domain. The computational domain is connected at the boundaries to two linear and homogeneous semi-infinite domains to accommodate the infinite extent of the railway track. The equations of motion for the three domains read where primes denote partial derivatives with respect to x, overdots denote partial derivatives with respect to t, and all parameters have been scaled by the beam's bending stiffness E I ; ρ 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 for the right semi-infinite domain; F k,W (x, w) and c d (x) are the scaled force exerted by the foundation and the scaled foundation damping of the computational domain, respectively; F 0 and v represent the scaled magnitude and the velocity of the moving load; δ(. . .) denotes the Dirac delta function, and H (. . .) the Heaviside function; 0 and L represent the positions of the left and right boundaries of the computational domain, respectively. Finally, w l , w r and w are the unknown displacements of the left and right semi-infinite domains, and of the computa- tional domain, respectively. The space and time dependency of the unknown displacements is omitted from the expressions for brevity.
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 as x tends to negative and positive infinity, respectively: lim For computational efficiency, the simulation is performed for the time interval when the moving load is close to and inside the transition zone. Therefore, the choice of initial conditions is crucial for ensuring that the infinite extent of the model is respected. In reality, before reaching a transition zone, the displacement field caused by a train can be considered to be in the steady state. Consequently, the initial conditions have to be chosen such that they represent the steady state at the beginning of the simulation. These initial conditions are referred to (in this paper) as the input initial conditions. The input initial conditions read where superscript "in" stands for the input initial conditions, which are derived in Sect. 2.3.1. Equations (1) to (11) constitute a complete description of the problem. In the next section, the solution method is presented.

Locally inhomogeneous and nonlinear system
Several time-domain methods are available for obtaining the solution to a system representing a railway track with nonlinear behaviour of the foundation [e.g., 28,37,39,40]. These methods are suitable for systems that exhibit nonlinear behaviour continuously throughout the simulation. An alternative method is using the Laplace transform sequentially, as shown by Hoving and Metrikine [41]. The main condition for this method to be applicable is that the system's behaviour is piecewise linear, implying that the system behaves linearly between the moments at which its parameters, being functions of the field variables (displacements, velocities, etc.), change abruptly (i.e., nonlinear events). This method has the potential of being computationally efficient for systems that have a limited number of nonlinear events.
For simplicity of the derivation in this section, the foundation's constitutive law is assumed to be bilinearly elastic. A more realistic constitutive relation, discussed in Sect. 3, is adopted for the results presented in this paper. Here, the force provided by the Winkler foundation springs is given by where k A d (x) and k B d (x) represent the foundation stiffness related to the first and second branches of the bilinear constitutive law (see Fig. 6, considering just branches A and B for loading and unloading), is the stiffness difference between the two branches, and w el represents the elastic displacement limit at which the stiffness changes from branch A to branch B. Note that, due to the inhomogeneity, both stiffness parameters k A d (x) and k B d (x) are functions of space; however, the elastic displacement limit w el is independent of the spatial coordinate. The choice of the parameters is discussed in Sect. 3.
Assuming that the system is in the linear regime at the start of the simulation and applying the Laplace transform over time to Eq. (2), the Laplace-domain equation of motion valid in the computational domain readŝ whereŵ 1 represents the unknown displacement in the Laplace domain; s = σ + iω is the Laplace variable, where σ is a small and positive real number and ω represents the angular frequency; subscript 1 represents that the analysis is performed for the system before the elastic displacement limit w el is exceeded for the first time. Furthermore,F ML 1 represents the Laplacedomain force exerted by the moving load andF IC 1 represents the forcing induced by the input initial conditions: It must be noted that the Dirac delta function which describes the moving load in the time domain becomes a spatially distributed load in the Laplace domain, as can be seen in Eq. 14.
Because of the spatial variation of the foundation stiffness and damping, Eq. (13) cannot be solved analytically for all stiffness and damping profiles. Therefore, the fourth-order spatial derivative is approximated using the Finite Difference Method. A central difference scheme of order O( x 6 ) is used inside the domain, and a hybrid between a central and a forward/backward scheme is used at the boundaries. The coefficients used for the finite difference discretization are given in Appendix A. As boundary conditions for the computational domain, the reaction forces delivered by the semi-infinite domains, namely the bending moment Eq. (6) and the shear force Eq. (7), are employed. These reaction forces are derived in Sect. 2.3.2 by imposing the displacement Eq. (4) and slope Eq. (5) of the computational domain as boundary conditions for the semi-infinite domains. Note that for the non-reflective boundaries derived in Sect. 2.3.2 to be correct, the computational domain must behave linearly at the boundaries. Therefore, the length of the computational domain must be chosen such that the expected nonlinear behaviour is located between its boundaries. After applying the Finite Difference Method to discretize the computational domain, the Laplace-domain equation of motion reads where K i j represents the bending stiffness matrix of the beam incorporating the contribution of the boundary conditions which is proportional to the unknown displacement, whileF D 1,i represents the contribution of the boundary conditions independent of the unknown displacement, which is regarded as an external forcing (see Sect. 2.3.2); I i j is the identity matrix. The Laplace-domain displacementŵ 1, j is obtained by leftmultiplying Eq. (16) by the inverse of the dynamic stiffness matrix (i.e., the term in the square brackets). Then, the inverse Laplace transform is numerically evaluated to obtain the solution in the time domain. Making use of the symmetry properties of the imaginary and real parts of the Laplace-domain spectrum, only positive frequencies are considered. For the results presented in this paper, the Trapezoidal Rule is used to evaluate the integral numerically.
The obtained time-domain solution is correct until the first nonlinear event, defined as the moment in time at which the solution exceeds the elastic limit w el at a certain location inside the computational domain. To obtain the correct solution after the first nonlinear event (i.e., time moment τ 1 ), the equation of motion of the computational domain, Eq. (2), is changed as follows: -The Winkler stiffness profile is modified by assigning the adequate stiffness to the nodes where the elastic limit has been exceeded. -A new time variable is introduced, namely t 2 = t − τ 1 for t ≥ τ 1 . Note that time t represents the global time. Furthermore, for clarity, the time moment of a nonlinear event in the global time axis t is represented with an overbar (τ n ), while in the local time axes t n , this moment is indicated without the overbar (τ n ). Figure 2 offers a graphical representation of the different time axes and nonlinear events. -To ensure continuity, the displacement and velocity of the previous system at τ 1 are prescribed as initial conditions for the new system: Note that these initial conditions differ from the input initial conditions defined in Sect. 2.3.1. -The boundary conditions are updated, as shown in Sect. 2.3.2. -The position of the moving load is updated accordingly.
The updated system behaves again linearly until the next nonlinear event. Therefore, the forward Laplace transform is applied with respect to the new time variable t 2 . The Laplace-domain equation of motion of the new system reads where s 2 represents the Laplace variable associated with the new time variable t 2 ,ŵ 2 is the unknown Laplace-domain displacement of the new system, and subscript 2 represents that the analysis is performed for the system in the second time interval. The termF ML 2,i is associated with the moving load acting on the new system, while the initial conditions for the new time interval are accounted for throughF IC 2,i : The Laplace-domain force exerted by the foundation is split into its contribution proportional to the unknown displacement, k d,2,iŵ2 , and the one independent of the unknown displacement, which is accounted for through the external forceF NL 2 , with superscript NL standing for nonlinear: The solution of the new system is obtained in the same manner as for the previous one. Moreover, the procedure of searching for the next nonlinear event, modifying the system and then solving it using the Laplace transform, is repeated. To this end, the procedure is generalized. The Laplace-domain equation of motion for the nth time interval reads where the generalized moving load and initial condition forces are given aŝ The described procedure is repeated until the whole solution is obtained in the time domain. The discretized displacement of the computational domain thus becomes where t is the time spacing, N is the index of the last time interval, and t max is the final moment in time of the simulation.
To obtain correct results, the continuity of displacements and velocities at the nonlinear events is of crucial importance. However, the Laplace-domain spectra of the two quantities exhibit a poor decay due to the applied initial conditions. Consequently, the numerical integration must be performed up to very high frequencies leading to a significant computational effort. In the following section, a method of incorporating the high frequencies without increasing the computational effort is presented.

Improvement of the frequency-spectra decay
The high-frequency regime of the Laplace-domain displacement, obtained from Eq. (24), is dominated by the initial conditions as follows: Similarly, the Laplace-domain velocityv n, j in the high- The initial displacement, velocity and acceleration in the numerators of Eqs. (28) and (29) are clearly independent of frequency. Therefore, the Laplace-domain spectra are dominated by the expressions in the denominator, exhibiting the slow decay proportional to 1/ω n . To incorporate the high frequencies without the need to integrate numerically up to very high frequencies, the high-frequency approximations, Eqs. (28) and (29), could be subtracted fromŵ n, j andv n, j , respectively. However, in doing so, a high peak is introduced in the remaining frequency spectra close to ω n = 0, which would require a very small step in frequency for obtaining accurate results. To overcome this issue, the highfrequency approximations are only subtracted over part of the frequency axis ω n ≥ ω A . The only criterion for choosing ω A is that it is sufficiently distant from the origin ω n = 0 to ensure that the resulting spectrum does not exhibit a high peak. The Laplace-domain displacement and velocity can now be expressed aŝ whereŵ imp n, j andv imp n, j represent the improved (i.e., with strong decay) Laplace-domain expressions of the displacement and velocity, respectively. To obtain the time-domain response, the inverse Laplace transform is evaluated numerically for the improved expressions, and analytically for the high-frequency approximations, as follows: where ω max represents the maximum frequency of integration, while I a (t n ) and I b (t n ) represent the timedomain images of the high-frequency approximations divided by the corresponding initial condition terms. Their expressions are not presented here for brevity; however, they can be easily computed using a symbolic mathematics tool (e.g., Maple). By evaluating the inverse Laplace transform analytically for the high-frequency approximations, frequencies up to infinity are actually included, which represents an improvement not just in computational time, but also in accuracy (although the results are still approximations sinceŵ imp n, j andv imp n, j are not completely zero at ω max ). Using the improvement presented in this section, the solution obeys the continuity conditions (i.e., displacement and velocity) imposed at nonlinear events. The input initial conditions and the nonreflective boundary conditions are derived in the next section.

Moving-load entry of the computational domain
At the start of the computation, the position of the moving load is at x = 0. For computational efficiency, this point must be as close as possible to the transition zone. To correctly represent the behaviour as described by Eqs. (1) to (3), which assume that the load has been active for a long time, the input initial conditions must be chosen such that the system is in the steady-state regime. Assuming that the system is initially in the linear regime (w(x, t = 0) < w el ), the input initial conditions are based on the steadystate field of the approaching load, also referred to as the eigenfield w e (x, t). The steady-state response of an infinite Euler-Bernoulli beam resting on a homogeneous and linear Winkler foundation to a moving constant load was derived using different methods (e.g., time-domain method, transform method) in, for example, [42]. Accounting also for the viscous damping of the supporting structure, the eigenfield reads where A 1 , A 2 , B 1 and B 2 represent the 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 B. It must be noted that w e (x, t) in Eq. (34) is real-valued. The eigenfield is derived assuming a homogeneous Winkler foundation. Consequently, choosing the input initial conditions based on the eigenfield is correct only if the field is not disturbed by the inhomogeneity (transition zone). The length of the computational domain and the position of the transition zone should thus be chosen such that the input initial state (i.e., displacement field, velocity field) based on the eigenfield has decayed before the inhomogeneity. Therefore, the part of the computational domain to the left of the transition zone has as input initial conditions the eigenfield part in front of the load, while the left semi-infinite domain has as input initial conditions the eigenfield part behind the load. Moreover, the input initial conditions at the inhomogeneity and to the right of it, including the right semi-infinite domain, are trivial. Equating t to 0 in Eq. (34) and in its time derivative, the input initial conditions, Eqs. (9) to (11), become With these input initial conditions, the system reaches the steady-state regime instantly at the start of the computation. Next, boundary conditions are derived such that the waves generated inside the computational domain are not reflected at the boundaries.

Derivation of the non-reflective boundary conditions
In Sect. 2.2, the boundary conditions for the computational domain, imposed bending moment Eq. (6) and shear force Eq. (7), were kept general. In this section, the reaction forces of the semi-infinite domains at the interfaces with the computational domain are derived. When these forces are prescribed as boundary conditions of the computational domain, the finite system will behave exactly as the infinite one. Therefore, these interface reaction forces constitute non-reflective boundary conditions for the computational domain. 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 t n to the equations of motion of the two semi-infinite domains, Eqs.
(1) and (3): whereŵ l,n andŵ r,n represent the unknown Laplacedomain displacements of the left and right semi-infinite domains, respectively, for the nth time interval, and F ML r,n represents the Laplace-domain moving load acting on the right domain, given by Eq. 26, but with a continuous spatial coordinate x. k l and k r represent the wavenumbers of the two semi-infinite domains and read where the branches of the complex wavenumbers are chosen such that Im(k h ) < 0 and Re(k h ) > 0. Furthermore,F IC l,n andF IC r,n represent the Laplace-domain initial condition forces given bŷ At infinity, the condition of zero displacements is imposed, Eq. (8), while at the interfaces the unknown Laplace-domain displacement and slope of the computational domain are prescribed: The Laplace-domain displacement of the two semiinfinite domains can be obtained by solving Eqs. (38) and (39) with the above-discussed boundary conditions. By taking the second-and third-order derivatives with respect to space and evaluating them at the interfaces, the reaction forces of the two semi-infinite domains are expressed as functions of the displacement and slope of the computational domain prescribed at the corresponding interfaces: where the entries of the matrices represent the dynamic stiffness coefficients giving rise to the boundary forces proportional to the unknown displacement and slope at the boundary; subscript V stands for shear force, M for bending moment, υ for translation, and φ for rotation. The coefficients read In addition, D l,n and D r,n are vectors containing the influence of the initial conditions Eq. (41) on the reaction forces, giving rise to boundary forces independent of the unknown displacement and slope of the computational domain. The two vectors are given by the following expressions: D r,n = ŵ r,n,p (L , s n ) w r,n,p (L , s n ) + k r,Vυ k r,Vφ k r,Mυ k r,Mφ ŵ r,n,p (L , s n ) w r,n,p (L , s n ) , whereŵ l,n,p (0, s n ) andŵ r,n,p (L , s n ) are the particular solutions that account for the initial condition forcing in Eqs. (38) and (39). Furthermore,V L,n andM L,n are the shear force and bending moment, respectively, exerted by the moving load on the right boundary after it has entered the right semi-infinite domain: Note that D r,n ,V L,n andM L,n as given in Equations (49) to (51), respectively, are valid for τ n−1 ≤ L v , meaning that the nonlinear events occur while the moving load is inside the computational domain. They still correctly describe the dynamics when the load is in the right domain, up to the moment a nonlinear event occurs. When nonlinear events occur while the moving load is in the right domain, this is divided into two domains, one behind the load and one in front, rendering the expressions for D r,n ,V L,n andM L,n lengthy for τ n−1 > L v . Therefore, these expressions are given in Appendix C.
To obtain the particular solutions in Eqs. (48) and (49), the Green's-function approach is used as follows:  48) and (49). Therefore, excitation variable ξ is smaller than or equal to the observation point x = 0 for the left domain and ξ is larger than or equal to x = L for the right domain. Consequently, the Green's functions are given by [43] Now, the only unknowns left for deriving the nonreflective boundary conditions are the initial conditions of the two semi-infinite domains in Eq. (41). The state (i.e., displacement field, velocity field) of the two domains at time moment τ n−1 consists of a superposition of the eigenfield and the waves generated inside the computational domain which have propagated to w f '(0,t) After computing the free-field slopes numerically, the boundary value problems needed to determine the state of the left and right domains (Fig. 3) are solved using the Laplace transform over global time t. Note that w f (0, t > τ n−1 ) and w f (L , t > τ n−1 ) are still unknown and by default equal to zero. Consequently, the discontinuity in the free-field displacements and slopes at τ n−1 introduces high-frequency content in its Laplace-domain counterparts.
To avoid this, an artificial smooth continuation is imposed on the free-field displacements and slopes for t > τ n−1 , which does not affect the response for t ≤ τ n−1 . Then, the Laplace-domain free-field displacements are given as follows: where C l , D l , C r and D r represent complex-valued amplitudes which read Note that the forward Laplace transform is applied with respect to the time variable t because the displacement and slope imposed as boundary conditions act from the time moment t = 0 until t = τ n−1 . Therefore,ŵ f (0, s) andŵ f (L , s) represent Laplace-domain history contributions for the new system and need to be computed for each nonlinear event.
To obtain the time-domain displacement and velocity of the two semi-infinite domains needed in Eq. (41) and thus for the derivation of the non-reflective boundary conditions, the inverse Laplace transform is applied to Eqs. (58), (59) and the corresponding velocities and is evaluated at τ n−1 : The particular solutions are now obtained by substituting Eqs. (64) to (67) in Eqs. (41), and (41) in Eqs.

End of the simulation
whereŵ e l,n,p (x, s n ) andŵ e r,n,p (x, s n ) represent the particular solutions accounting for the parts of the eigenfield in the left and right domains, respectively, at τ n−1 . The integration over ξ can be performed analytically, while the inverse Laplace transform evaluated at τ n−1 should be performed numerically. Moreover, the spatial derivatives needed in Eqs.
The contribution of the boundary conditions which is proportional to the yet unknown displacement and slope is incorporated into the beam's bending stiffness matrix K i j (see Eq. (24)), while the contribution which is independent of the unknown displacement and slope is accounted for through the boundary-forcing vectorF D n,i . As can be seen, the beam's bending stiffness matrix does not change from one system to the other; however, the boundary-forcing vector needs to be updated at each system change. A synoptic chart of the developed algorithm is presented in Fig. 4.
Although the obtained non-reflective boundary conditions have been derived for the moving load problem, the same procedure can be applied for any arbitrary loading. Next, the Winkler foundation model used for the results presented in this paper is described.

Winkler foundation
A variety of models have been used to represent the supporting structure in 1-D systems for railway applications (e.g., linear elastic, bilinear elastic, cubic superlinear, etc.); a good overview of models for the foundations in 1-D systems can be found in [44]. One of the most used nonlinear elastic models is the cubic superlinear one [e.g., [30][31][32][33]35]. It assumes that the foundation springs exert a reaction force (visualized in Fig. 5) proportional to the displacement, through a linear stiffness term k l , and one proportional to the displacement cubed, through a nonlinear stiffness term k nl . In this paper, the constitutive relation of the Winkler foundation is based on the cubic super-linear model, but also incorporates the possibility of plastic deformation by selecting a different unloading path than that of the loading, as seen in Fig. 6.
The loading path of the chosen constitutive relation approximates the cubic super-linear model through a piecewise-linear profile (Fig. 5) to accommodate the solution method presented in Sect. 2.2. Firstly, k A d is assumed to be equal to the stiffness of the equivalent linear model [45]. This assumption is based on the ballast being relatively well compacted at the start of the simulation, represented in Fig. 5 through the nonzero displacement at zero force in the piecewise-linear approx-   (2) imation. This implies that the initial soft response of the foundation, as sometimes encountered, is excluded. Furthermore, assuming that the compaction is uniform along the track, the compacted configuration is taken as reference (zero displacement for zero force in Fig. 6). Secondly, the steady-state eigenfield is assumed to be in the linear regime. This assumption is based on the fact that in the homogeneous parts of the railway track, the steady-state displacement field induced by a train does not lead to significant degradation of the supporting structure. Thirdly, the elastic displacement limit w el is chosen relative to the eigenfield's maximum displacement w e max in the soft part of the computational domain, where the ratio q = w el /w e max is larger than 1. If w el is not exceeded during the simulation, the system remains in the linear regime (branch (1) in Fig. 6). At locations where w el is exceeded, the corresponding part of the foundation enters the second loading branch k B d . The value of k B d is chosen such that it approximates the cubic super-linear model (Fig. 5).
The parameters of the cubic super-linear model which is approximated are chosen as similar to the ones used in other publications [i.e., 31,32]. However, the parameters for the unloading path, k C d and k D d , are not well known. For the moment, the parameters are chosen such that the overall constitutive relation resembles the results of cyclic loading experiments on granular material [i.e., 46], but specific additional experiments or 3-D modelling might be needed in order to choose these parameters realistically. The model developed in this paper can even incorporate cyclic behaviour, and the parameters of the constitutive relation can be modified for each cycle to accommodate the behaviour as observed in the literature: changing stiffness with increasing cycle number, as well as a changing area (energy dissipation). However, the computational time required for the simulation of thousands of trains passages will be high. Nonetheless, valuable observations can be made from the short-term cyclic behaviour which have relevant implications for the long-term behaviour.
The constitutive model incorporates the possibility of separation between the rail and the supporting structure (Fig. 6). When the displacement of the beam is larger than the plastic deformation (w > w pl ), the separation of the beam and foundation occurs. This results in a zero foundation force F k,W as depicted in Fig. 6. However, this is only allowed at the location where the plastic deformation has been activated; in the parts of the computational domain with no plastic deformation, the beam is in permanent contact with the supporting structure. In case of separation, besides the foundation stiffness, also the foundation damping is modified. Consequently, the foundation damping is given by where w pl represents the plastic deformation which has a negative value and c S d is the remaining damping coefficient in case of separation which, in fact, could represent part of the internal damping of the beam.
The parameters of the constitutive relation are not only functions of the displacement, but also functions of space, due to the inhomogeneity of the modelled system. To study the influence of the transition smoothness on the plastic deformation in the transition zone and on the radiated wave field, the spatial profile of the foundation stiffness is chosen based on a sine squared as follows: where h = {A, B, C, D}; x tc represents the transition centre, l t the transition length, and p 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 input initial conditions, as discussed in Sect. 2.3.1. In addition, a similar spatial variation is chosen for the foundation damping. In the rest of the paper, the damping is expressed through the ratio ζ which is defined similar to that of a singledegree-of-freedom system: Therefore, by maintaining a constant damping ratio ζ throughout the system, the spatial variation of the foundation damping is proportional to the square root of that of the stiffness (branch A), except for the parts of the beam which have lost contact with the foundation.

Results and discussion
Here, the proposed model is first validated by considering a limit case and comparing the obtained results to a semi-analytical solution. Then, the time-domain displacement field is presented for two specific cases and the influence of the nonlinear foundation on the transition radiation is highlighted. Afterwards, the influence of the transition length, load velocity and stiffness dissimilarity on the plastic deformation is assessed through a parametric study. Finally, the influence of the nonlinear foundation on the radiated energy and on the energy input is discussed. The parameters which are kept constant throughout the presented results are given in Table 1, while the ones which are varied are mentioned for each case individually.

Validation and convergence
To validate the solution derived in Sect. 2, a limit case is considered, in which the foundation is piecewise homogeneous and behaves linearly, but for which artificial nonlinear events are introduced in the solution. To this end, a soft-to-stiff case is considered where the foundation stiffness coefficients k  Stiff-soft stiffness ratio p 5 Elastic displacement limit ratio q 1. The limit-case solution is compared to the semianalytical solution of a piecewise-homogeneous and linear system. This semi-analytical benchmark solution can be obtained by solving a system of two semi-infinite domains with different foundation stiffnesses by using the Fourier transform over time. In the Fourier domain, after imposing interface conditions and the condition of zero displacements at infinity, the displacements of the two domains can be obtained analytically [i.e., 16,42]. To obtain the time-domain solution, the inverse Fourier integral can be applied numerically.
The error e(x) presented in Fig. 7 is defined as the summed-over-time absolute value of the difference between the limit-case solution w lin and the benchmark solution w bench , divided by the summation of the absolute value of the benchmark solution over time: This error is caused by two main factors: the sequential application of the Laplace transform and the finite difference discretization. The left panel of Fig. 7 presents the total error. To isolate the error caused by the finite difference discretization, the case with no artificial nonlinear events is presented in the right panel of Fig. 7.
To test the convergence of the derived solution, the maximum frequency has been varied (also done in the numerical integration for the benchmark solution). Note that by changing the maximum frequency, according to the Nyquist sampling rule, the time stepping also changes. Figure 7 shows that the solution derived in Sect. 2 converges to the correct one as the maximum frequency increases. The higher relative error in the stiffer part of the computational domain can be explained by the smaller displacements. A higher maximum frequency leads to smaller error; however, the computational effort increases significantly. For the rest of the results presented in this section, the maximum frequency was chosen as 1000 Hz.

Displacement field in the time domain
To study the effect of the nonlinear foundation on the wave field excited during the transition radiation process, a relatively severe case is presented. The load velocity is chosen as 95% of the minimum phase velocity in the soft part of the computational domain (c l,ph,min ), and the transition length l t is chosen as 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 35 m. The influence of the nonlinearity is highlighted by comparing the response to the linear case, as shown in Fig. 8. (d) t = 0.097 s Fig. 9 Time-domain displacement field for the nonlinear system (solid line), the plastic deformation (dash-dot line) and the elastic displacement limit (light grey; given in Table 1)  The plastic deformation close to the transition zone is a result of constructive interference of the approaching eigenfield and the waves generated in the transition zone. As in the linear case, the eigenfield interacts with the transition zone, generating waves which mostly propagate towards the softer medium. By constructively interfering, the displacement under the load exceeds the elastic displacement limit (indicated by the grey line in Fig. 8) giving rise to permanent deformation in the foundation to the left of the transition zone.
In the nonlinear system, the free field has both larger amplitude and is sustained for a much longer period of time when compared to the linear one. This can clearly be seen in Fig. 8. Both the larger amplitude and longer duration are consequences of the beam's separation from the foundation. The loss of contact leads to larger upward displacement which in turn affects the wave field in areas without loss of contact. The loss of contact also causes the loss of external damping which results in the longer duration of the vibration. Moreover, the shape of the radiated waves is also changed. In the nonlinear case, the contact loss leads to a slight increase in the wavelength of the free field ( Fig. 8 panels  (d), (e) and (f)). This is in accordance with the findings in Sect. 4.4.
When the moving load travels from a stiff medium to a soft one, the displacement field does not exceed the elastic displacement limit (at all), meaning that plastic deformation does not develop, as seen in Fig. 9. As in the previous case (Fig. 8), the large-amplitude free waves propagate into the softer medium, which in this case is in the forward direction. This leads to much less constructive interference between the eigenfield and the free field and thus to no plastic deformation. Therefore, in the remainder of the paper, only the softto-stiff case is presented. However, if the load travelled super-critically (in the soft medium or in both media) for the stiff-to-soft case, it would move faster than the minimum phase velocity of the free waves, probably leading to more pronounced constructive interference, which could induce plastic deformation to the right of the transition zone.

Parametric study
In this section, the damage occurring in the supporting structure is addressed as a function of the load velocity, the transition length and the stiffness dissimilarity p through a parametric study. The area of the plastic deformation, A pl = w pl (x) dx, is chosen as the quantifier for the damage in the supporting structure. The three parameters chosen to be varied influence the transition radiation phenomenon most. Furthermore, these parameters can be adjusted in the design stage of railway tracks to minimize damage in the supporting structure. The plastic deformation area versus varying transition length is presented in Fig. 10 for different load velocities and damping ratios. The plastic deformation area decreases as the transition length increases, as could be expected, and the decreasing trend is quasimonotonic for all load velocities and damping ratios. For the damping ratio ζ = 0.05 (left panel in Fig. 10), the transition-length range in which plastic deformation occurs increases with increasing velocity until the minimum phase velocity, beyond which it decreases. However, for the super-critical case (v = 1.05 c l,ph,min ) and small transition lengths (l t = 0, . . . , 3 m), the plastic deformation is still larger than for the critical case (v = 1.00 c l,ph,min ), which makes the two lines intersect. This larger plastic deformation is caused by the fact that parts of the foundation experience additional plastic deformation produced in the second loadingunloading cycle. This is not the case for larger transition lengths (l t = 4, . . . , 6 m), and because the displacement of the eigenfield under the load is smaller in the super-critical case than in the critical case, a smaller plastic deformation area results, which explains the intersection of the two lines. Furthermore, for ζ = 0.20 (right panel in Fig. 10), additional plastic deformation is not produced in the second loading-unloading cycle for any of the transition lengths due to the higher damping ratio, and the range as well as the plastic deformation area just decrease as the velocity increases beyond the minimum phase velocity.
Furthermore, the analysis shows that it is not just the duration of passage t p = l t v that governs the resulting plastic deformation (which could be intuited), but also the absolute values of the transition length l t and load velocity v separately. It can clearly be seen in Fig. 10 that for the same duration of passage, significantly different plastic deformation values are observed for different load velocities.
In Fig. 11, the plastic deformation area is presented as a function of the velocity of the moving load for different transition lengths and damping ratios. The plastic deformation area increases with increasing velocity until close to the minimum phase velocity, beyond which it decreases. For ζ = 0.05 (left panel in Fig. 11), the critical velocity appears to be around 1.05 c l,ph,min for l t = 0.1 m, and its value decreases with increasing transition length reaching 1.00 c l,ph,min for l t = 4 m. This shows that the critical velocity for the plastic deformation area is dependent on the transition length. Moreover, the increasing trend (sub-critical velocities) and the decreasing trend (super-critical velocities) have different slopes, which is explained by the fact that in the super-critical cases parts of the foundation experience additional plastic deformation caused in the second loading-unloading cycle. For ζ = 0.20 (right panel in Fig. 11), a similar behaviour to the case of ζ = 0.05 is observed. However, the magnitude of the plastic deformation area is smaller and the slopes of the increasing and decreasing trends (left and right of the critical velocity) are almost identical because the foundation does not experience additional plastic deformation caused in the second loading-unloading cycle.
The plastic deformation area as a function of the stiffness dissimilarity p is presented in Fig. 12 for different load velocities and transition lengths. The plastic deformation area increases with increasing stiffness dissimilarity for all velocities and transition lengths presented, as could be expected. The increasing trends tend to constant values, which could be obtained in It is clear from the trends that the maximum sensitivity of the plastic deformation area occurs at small stiffness dissimilarities. Therefore, the maximum gain in decreasing the damage of the supporting structure can be obtained at small stiffness dissimilarities. For small stiffness dissimilarity, the maximum plastic deformation area is observed for a load velocity v = 1.00 c l,ph,min , but for larger stiffness dissimilarity the maximum plastic deformation area is obtained for a load velocity v = 1.05 c l,ph,min . This shows that the critical velocity, when it comes to the plastic deformation area, is dependent also on the stiffness dissimilarity (next to transition length).
To conclude, studying the influence of the transition length, load velocity and stiffness dissimilarity on the plastic deformation area provides valuable information about the value ranges of these parameters where the initial design of transition zones should aim at so as to minimize the damage in the supporting structure. Next, the transition radiation phenomenon is studied from an energy point of view.

Energy radiation
In this section, the influence of the nonlinearity on the transition radiation is studied from the energy point of view, which gives additional insight into the properties of the radiated field. To investigate the energy radiation solely due to the transition radiation phenomenon, the study is restricted to sub-critical load velocities. Based on the solution derived in Sect. 2.2, the energy flux through cross sections of the beam to the left and right of the transition zone can be computed as seen in [10,16,42]. To study the energy associated with the radiated waves induced by the load passing the transition zone, the power flux corresponding with the free field S f (x, t) through a certain cross-section x of the Euler-Bernoulli beam is considered [42]: where w f = w − w e represents the free-field displacement; w is the displacement given by Eq. 27, and w e is the eigenfield displacement given by Eq. 34. It must be noted that in Eq. (76) the power flux propagating in positive x-direction is considered; to obtain the power flux in negative x-direction, a minus sign needs to be added. The total free-field energy flux E f (x) through a cross section is found by integrating the power flux over time. Moreover, to visualize the spectral energy density corresponding to the free field, E f (x) can be rewritten in terms of the Fourier-domain counterparts of the time-domain quantities, as shown in [15,16]: where the integrand in the last expression represents the spectral energy density E s (x, ω), ω represents the Fourier-domain angular frequency (rather than that in the Laplace domain as introduced in Sect. 2.2), the tilde represents the Fourier-domain quantities,ṽ represents the Fourier-domain velocity, and the asterisk represents complex conjugation. It must be noted that Eq. (77) does not capture the entire transition radiation energy; it only covers the energy carried by the isolated free field. The total energy flux through a certain cross section contains the eigenfield contribution, the free-field contribution, as well as an interference term; the last two both contribute to the radiation energy, while only the free-field contribution is considered here. Furthermore, the computed energy does not represent the total free-field energy because part of it is absorbed by the foundation before reaching the considered cross sections left and right of the transition. Moreover, due to the spatial variation of the foundation damping, the energy propagating in the stiff domain is damped more than that propagating in the soft domain. To avoid the last issue, for the computations performed in this section, the spatial damping profile is maintained constant throughout the computational domain (where ζ is defined with respect to the soft domain), except for the parts where the beamfoundation separation occurs Eq. (72).
The influence of the transition length on the freefield energy is now addressed through a parametric study. The energy associated with the leftwardpropagating free field E f left , presented in the left panels of Fig. 13, decreases with increasing transition length, as could be intuited. However, for small transition lengths, E f left in the lightly damped case (panel a in Fig. 13) is smaller in the nonlinear system as compared to that in the linear one. This can be explained by the fact that part of the energy is consumed to plastically deform the foundation, but also by stronger radiation in the rightward direction (panel b in Fig. 13). Furthermore, although E f left decreases with increasing transition length, the decrease rate is smaller in the nonlinear case as compared to that in the linear one. Consequently, for some values of the transition length (l t ≈ 3, . . . , 7 m in panel a of Fig. 13), E f left is larger in the nonlinear case. However, the behaviour described above is not general. The bottom-left panel in Fig. 13 shows that in the system with higher damping ratio (ζ = 0.20), E f left is higher in the nonlinear system than in the linear one for all transition lengths.
Furthermore, the energy associated with the rightward-propagating free field E f right (presented in the right panels of Fig. 13) is significantly higher in the nonlinear case as compared to that in the linear one.
Moreover, E f right is not largest for the shortest transition length, but for a transition length l t ≈ 0.5, . . . , 2.5 m (depends on the load velocity). This is specific to the nonlinear case since for all the linear cases considered, the free-field energy decreases with increasing transition length. In addition, for the considered cases, E f right is clearly smaller than E f left , implying that the free field radiated into the soft part of the system carries most energy, although this is not a necessity as shown in [16].
The spectral energy density E s for the same system as in Fig. 8  shift towards the lower frequencies as well as a decrease in energy in the nonlinear case, when compared to the linear system. The lower frequencies of the radiated waves in the nonlinear system can also be observed in the time-domain response (Fig. 8, panels d), e) and f)) through the larger wavelengths of the left-propagating waves. Furthermore, the energy propagating rightward E s right in the nonlinear system has, next to a higher magnitude (as already observed in the right panel of Fig. 14), also higher-frequency content as compared to the linear system.

Energy input
Besides the influence of the nonlinear foundation on the energy radiation, presented in the previous subsection, its influence on the energy input from the moving  load could also offer valuable insight. While the energy radiation refers to the far field, the energy input offers information about the near field, as it relates to the contact point between the moving load and the beam. More energy input leads to more energy radiated and/or more energy dissipated into the foundation, the latter leading to damage. Therefore, the energy input or the maximum power input could represent a good indicator of the potential damage occurring in the foundation. To investigate the energy input solely due to the transition radiation phenomenon, the study is restricted to sub-critical load velocities.
The energy input is defined as the power input integrated over time. Due to the damping present in the structure, the power input P in is nonzero over the whole time axis (i.e., even when the load is not in the vicinity of the transition zone), and therefore, the energy input is infinite. Consequently, the difference in energy input between the linear and nonlinear cases is presented instead. The difference in energy input E in is given by the following expression: where P in lin andẇ lin represent the power input and velocity at the contact point, respectively, of the linear case. Figure 15 presents the difference in energy input E in as a function of the transition length for different velocities and damping ratios. It can be observed that the maximum difference in energy input does not occur at the smallest transition zone, reinforcing the findings in the energy radiation study (right panels in Fig. 13). Furthermore, the difference in energy input increases with increasing velocity, but the difference is smaller in the higher-damping case (right panel in Fig. 15).
The normalized maximum power input P in max is presented in Fig. 16 as a function of the transition length for different velocities and damping ratios. The maximum power input is normalized by the steady-state power input in the soft domain P in soft . It can be observed that the maximum power input decreases with increasing transition length and decreasing velocity of the moving load. Moreover, the maximum power input in the transition zone can be a factor 1.8 larger than the one in the steady state. Part of this additional power input could cause damage in the supporting structure.
The difference between the linear and the nonlinear cases (for the first passage of the moving load) is small and cannot be seen in Fig. 16, even though the difference in total energy input can be significant (Fig. 15). However, the increase is much more drastic for the second passage of the moving load (i.e., of the system that has already been deformed plastically) as can be seen in Fig. 17. Both the maximum power input (left panel in Fig. 17) and the difference in energy input (right panel in Fig. 17) exhibit a considerable increase compared to the first load passage. The maximum power input and the energy input have the potential to be good indicators of the damage occurring in the foundation of the railway track. However, more extensive research into these two indicators needs to be performed in order to justify this.

Conclusions
In this paper, the influence of the foundation's nonlinear behaviour on the waves generated by a moving load crossing an inhomogeneity in the foundation as well as the resulting plastic deformation have been studied. To this end, a one-dimensional model has been formulated, consisting of an infinite Euler-Bernoulli beam resting on a locally inhomogeneous and nonlinear Winkler foundation, subjected to a moving load. The reaction of the Winkler foundation has been characterized by a piecewise-linear (in displacements) constitutive relation which accounts for perma-nent deformations. The foundation's piecewise-linear behaviour allows to obtain the solution by sequentially applying the Laplace transform over time, while the Finite Difference Method has been used for the spatial discretization. The infinite extent of the system has been accounted for through a set of non-reflective boundary conditions, derived by replacing the semiinfinite domains by their response at the interfaces, and through the input initial conditions based on the steady-state response of a beam with homogeneous foundation subject to the moving load. The solution has been validated for the limit case of a linear and piecewise-homogeneous foundation against a semianalytical benchmark solution.
Results have shown that the plastic deformation originates from constructive interference between the waves excited at the transition and the response to the load's deadweight. Furthermore, through a parametric study conducted for the soft-to-stiff transition, it has been found that the plastic deformation decreases quasi-monotonically with increasing transition length, and that the decrease rate depends on the velocity of the load and on the magnitude of the foundation damping. Moreover, for super-critical velocities, additional plastic deformation is generated in the second loadingunloading cycle of the foundation for small transition lengths and/or large stiffness dissimilarities. The critical velocity related to the plastic deformation area, which was chosen to quantify the damage in the supporting structure, was observed to be dependent on the transition length and on the stiffness dissimilarity. In addition, the resulting plastic deformation area has been observed to be influenced not only by the time of passage of the transition zone, defined as the transition zone length divided by the load's velocity, but also by these quantities individually. Finally, when the load travels from a stiff medium to a soft one, no plastic deformation has been observed, which is due to less pronounced constructive interference.
The influence of the foundation's nonlinear behaviour on the radiated field has also been studied from the energy point of view, by considering the energy flux associated with the free field through cross sections left and right of the transition zone. Results have shown that in the nonlinear system, the maximum rightward-propagating energy flux does not occur at the smallest transition length, but at a larger one, a finding which is reinforced by the difference in energy input between the linear and nonlinear case. This feature is specific to the nonlinear system since in the linear case, the energy flux has been observed to always decrease with increasing transition length. Furthermore, the spectral energy densities have not only shown that the nonlinearity redistributes the energy between frequencies, but have also highlighted the redistribution between the soft and stiff media. Moreover, it has been observed that the maximum power input and the difference in energy input between the linear and nonlinear case drastically increases for the second passage of the moving load (i.e., of the system that has already been deformed plastically). This suggests the use of these energy quantities as possible indicators of the damage in the supporting structure.
Although one-dimensional models are not able to model all phenomena, they are useful for initial assessments. The model presented here can be used for preliminary designs of transition zones in railway tracks. Given the stiffness dissimilarity and/or the magnitude of the initial plastic deformation, the optimum length of the transition zone and the maximum velocity of the train can be obtained such that the damage in the railway track is minimized.
where the nodes −1, 0, M + 1 and M + 2 represent the so-called ghost nodes. These nodes are not part of the computational domain; therefore, the displacements of these nodes (ŵ −1 ,ŵ 0 ,ŵ M+1 ,ŵ M+2 ) are expressed in terms of the displacements of the nodes inside the computational domain by using the four nonreflective boundary conditions. The resulting equations are included in Eq. (24).
x = x P , namely the continuity of displacement and