Transition radiation in a nonlinear and infinite one-dimensional structure: a comparison of solution methods

Transition zones in railway tracks are locations with a significant variation of track properties (i.e. foundation stiffness) encountered near structures such as bridges and tunnels. Due to strong amplification of the track’s response, transition zones are prone to rapid degradation. To investigate the degradation mechanisms in transition zones, researchers have developed a multitude of models, some of them being very complex. This study compares three solution methods, namely an integral-transform method, a time-domain method, and a hybrid method, with the goal of solving these systems efficiently. The methods are compared in terms of accuracy, computational efficiency, and feasibility of application to more complex systems. The model employed in this paper consists of an infinite, inhomogeneous, and piecewise-linear 1-D structure subjected to a moving constant load. Although the 1-D model is not particularly demanding computationally, it is used to make qualitative observations as to which method is most suitable for the 2-D and 3-D models, which could lead to significant gains. Results show that all three methods can reach similar accuracy levels, and in doing so, the time-domain method is most computationally efficient. The integral-transform method appears to be efficient in dealing with frequency-dependent parameters, while the time-domain and hybrid methods are efficient in dealing with a smooth nonlinearity. For multi-dimensional models, if nonlinearities and inhomogeneities are considered throughout the depth, the time-domain method is likely to be most efficient; however, if nonlinearities and inhomogeneities are limited to the surface layers, the integral-transform and hybrid methods have the potential to be more efficient than the time-domain one. Finally, although the 1-D model presented in this study is mainly used to assess the three methods, it can also be used for preliminary designs of transition zones in railway tracks.


Introduction
Transition radiation occurs when a source moves with a constant velocity through or in the vicinity of an inhomogeneous medium [1,2]. It occurs, for example, when a train crosses a transition zone, which in railway tracks is an area with substantial variation of track properties (e.g. foundation stiffness) encountered near rigid structures such as bridges, tunnels or culverts. In railway tracks, transition radiation is emitted in the form of elastic waves that can constructively interfere with the deformation field caused by the dead weight of the moving train. Consequently, transition radiation has been identified 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][7]. This leads to transition zones requiring 3-8 times more frequent maintenance than the regular parts of the railway track [8,9]. In order to improve the design of transition zones and propose mitigation measures, engineers and researchers have developed a multitude of models to describe and predict the response of such systems. Some researchers focus on modelling the ballast behaviour (e.g. [10][11][12][13][14][15]), which is of crucial importance for the degradation of transition zones in ballasted tracks, some focus on transition radiation itself, while others combine both [7,16,17]. This paper belongs in the last category.
The first study on transition radiation of elastic waves was published by Vesnitskii and Metrikin [18]. It considered an infinite string on a piecewisehomogeneous Winkler foundation subjected to a moving point force and a moving point mass. The response of the system was obtained using the method of images for both types of moving loads, while for the moving point mass the Fourier transform over space was additionally used. To also account for the flexural rigidity of the system, Vesnitskii and Metrikin [19] analysed transition radiation in a semi-infinite beam on elastic springs clamped at one end. Once again, the method of images combined with the Fourier transform over space was used to obtain the solution. Later on, the problem of a beam resting on a piecewise-homogeneous Winkler foundation and subjected to a moving load was solved using modal expansion techniques [20,21]. To study wave propagation in the ground caused by transition radiation, 2-D models of a piecewise-homogeneous continuum were analysed using the Fourier transform [22], and the Fourier transform combined with mode matching [23,24]. Lately, a combination of the Fourier transform over time and the Fourier expansion over space was used to study transition radiation in a discretely supported Timoshenko beam [25]. Also, transition radiation was studied in a piecewise-linear 1-D system using a sequential Laplace transform method combined with a finite difference discretisation of the spatial dimension [7]. Furthermore, the interaction between a moving oscillator and an inhomogeneous and infinite structure was analysed by means of the Green's function method [26], which can be considered as a hybrid between an integral-transform method and a time-domain method.
Although at the beginning researchers have mainly used integral-transform methods to study the dynamics of elastic structures subjected to moving loads, the development of computers has led to a shift towards numerical methods such as the finite element method combined with time-integration methods (e.g. the Newmark method [27][28][29]). One advantage of these timedomain methods is that the geometry of the transition zone can be modelled accurately [6,[30][31][32][33][34][35]. Another advantage is that the nonlinear behaviour of the foundation or the nonlinear interaction between the vehicle and the structure can easily be handled by timedomain methods [36][37][38][39][40]. One of the disadvantages of the standard time-domain methods is that the system must be finite to be solved numerically, while the railway track is practically infinite, potentially causing artificial reflections at the boundaries of the finite domain. To overcome this problem, absorbing boundaries (e.g. perfectly matched layers) have been used in numerous studies (e.g. [41]). Another challenge caused by the structure being finite is that the vehicle's action on the structure is of finite duration causing an unrealistic transient behaviour at the entrance and exit of the vehicle. This has been elegantly solved by analysing the system in the moving reference frame, approach which is sometimes called the moving element method (e.g. [42,43]).
As it can be seen, a multitude of methods have been applied by researchers and engineers to investigate transition radiation in railway applications. Most methods fall into three main categories, namely integraltransform methods, time-domain methods, and hybrid methods (i.e. a combination of integral-transform and time-domain methods). The advantages and disadvantages of these methods have partially been discussed for each method independently, but without much direct comparison between them. This paper aims at comparing three solution methods, one corresponding to each category, for analysing transition radiation in a 1-D model consisting of an infinite Euler-Bernoulli beam resting on an inhomogeneous and piecewiselinear Kelvin foundation subjected to a moving constant load. The integral-transform method, namely the sequential Laplace transform method [7], assumes that the system's behaviour is piecewise linear and deals with the linear parts of the solution in the Laplace domain. To accommodate the inhomogeneity, the finite difference method is used for the spatial discretisation, and non-reflective boundary conditions are imposed to ensure the infinite extent of the structure. The timedomain method is the more conventional approach; the solution is obtained by applying the finite element method for the spatial discretisation, combined with a set of non-reflective boundaries, and the Newmarkβ as the time-stepping method. The hybrid method, namely the pseudo-force method [44,45], treats the nonlinearity and part of the inhomogeneity as external forces resulting in a fictitious linear and piecewisehomogeneous base structure. Taking advantage of the linearity, the Green's function of the base structure is obtained (the integral-transform part of the method), and the response is expressed through convolution integrals of the Green's functions and the forces that account for the nonlinearity and part of the inhomogeneity (the time-domain part of the method). Finally, because these forces are state dependent, the relation for the response at each time step is implicit, and it is therefore solved iteratively.
The three solution methods are compared in terms of accuracy, computational efficiency, and feasibility of application to more complex systems. In this paper, the feasibility of application to more complex systems refers to the feasibility of the methods to deal with frequency-dependent properties of the structure, to deal with a smooth nonlinearity (as opposed to the piecewise-linear one), and to apply the solution methods to 2-D and 3-D models. It must be emphasised that although the 1-D model described is not particularly demanding from the computational point of view, the aim of the comparison is to establish the most suitable method to be applied to multi-dimensional models, where the choice of the proper solution method could lead to significant gains.
The novelty of the current study is threefold. Firstly, a thorough and direct comparison between integral-transform, time-domain, and hybrid methods to describe the behaviour of an infinite and nonlinear system has not been presented in the literature. Such a comparison enables engineers and researchers to choose the most suitable solution method for solving the specific problem they are facing. Secondly, the application of the pseudo-force method to analyse transition radiation in a nonlinear and infinite structure is presented here for the first time. Thirdly, the non-reflective boundary conditions formulated for the time-domain solution method, which enable the finite domain to behave exactly as the infinite one, are derived analytically; to the best of the authors' knowledge, it is for the first time this analytical approach is used to solve a moving load problem containing changes in the foundation properties (see Sect. 2.4 for more details).
Finally, it must be noted that although the 1-D model presented in this study is mainly used to assess the three methods, it can also be used for preliminary designs of transition zones in railway tracks.

Problem statement
The model formulated in this section is composed of a constant moving load acting on an infinite Euler-Bernoulli beam that rests on a smoothly inhomogeneous and nonlinear Kelvin foundation and is depicted in Fig. 1. The model is divided in three domains: the computational domain that is nonlinear and inhomogeneous, and two linear and homogeneous semi-infinite domains to accommodate the infinite extent of the railway track (Fig. 1). The combined equation of motion for the three domains reads x ≥ L , where∂ ∂ x and∂ ∂t represent generalised derivatives with respect to space and time, respectively. All parameters in Eq. (1) have been scaled by the bending stiffness E I of the beam; ρ is 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 semiinfinite domain, respectively; k d,r and c d,r represent the same quantities for the right semi-infinite domain; f k,c (x, w c ) and c d,c (x, w c ) represent the scaled force and x e ≤ 0 represents the position of the moving load at t = 0. Also, w l , w r and w c represent the 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. Furthermore, the ≤ and ≥ signs in the definition of w(x, t), f k (x, w), and c d (x, w) for all domains emphasise that there is continuity in these quantities at the interfaces between the three domains.
At the interfaces between the three domains, continuity in displacement and slope as well as in shear force and bending moment is imposed. Furthermore, due to the presence of damping, the displacements at infinite distance from the moving load are zero. This implies that the left and right domains are at rest as x tends to negative and positive infinity, respectively. The interface and boundary conditions thus read where primes denote classical partial derivatives with respect to x.
As for the initial conditions, unless x e is placed far away from the transition (which is undesirable because it leads to an increase in computational time), the choice of the initial state (i.e. displacement and velocity fields) affects the response in the computational domain. The response caused by a train before it reaches a transition zone can be considered to be in the steady state. Therefore, the initial conditions should be selected such that the response at the start of the simulation is in the steady-state regime, and are thus based on the steadystate field, also referred to as the eigenfield w e (x, t). The eigenfield generated by a moving constant load acting on an infinite Euler-Bernoulli beam that rests on a homogeneous Kelvin foundation reads [7,46] wherex = x − x e , A 1 , A 2 , B 1 , and B 2 are complexvalued amplitudes, and k e 1 , k e 2 , k e 3 , and k e 4 are the complex-valued wavenumbers of the eigenfield. Their expressions are given in Appendix A in [7]. It must be noted that w e (x, t) in Eq. (7) is real-valued and it refers to the eigenfield in the left domain throughout the paper.
Imposing the initial state based on the eigenfield is correct only if the field is not disturbed by the inhomogeneity (transition zone) because the eigenfield is derived by assuming a homogeneous foundation. The location x e of the moving load at t = 0 should thus be chosen such that the initial state based on the eigenfield has decayed before the inhomogeneity. Consequently, where overdots denote classical partial derivatives with respect to t. As for the nonlinear behaviour of the foundation, the constitutive relation from Ref. [7] is chosen, which is summarised in the following. This constitutive model approximates the observed response of a granular material to cyclic loading (e.g. [47]) through a piecewiselinear profile (Fig. 2). While the loading path of such a material is relatively well known (a multitude of studies use the cubic superlinear model), the parameters of the unloading path (i.e. k C d,c and k D d,c in Fig. 2) are less well known and are therefore chosen such that the overall constitutive relation resembles the response of granular materials to cyclic loading. For the loading path, to approximate the cubic superlinear model by a piecewise-linear profile, a few assumptions are made. Firstly, the initial branch of the loading path with stiff-ness k A d,c is chosen to be the same as the stiffness of the equivalent linear model [48]. This assumption is correct provided that the ballast has been well compacted. Secondly, the eigenfield is assumed to be in the initial branch with stiffness k A d,c . This assumption comes from the fact that accelerated degradation is not observed in the homogeneous parts of the railway track where the response is in the steady state. Finally, the elastic displacement limit w el (Fig. 2) is assumed to be a fraction of the eigenfield's maximum displacement w e max in the soft part of the track, and, according to the second assumption, the ratio q = w el /w e max must be larger than 1.
Path 1 in Fig. 2 represents the constitutive behaviour for locations in the computational domain where w el is not exceeded during the simulation. If w el is exceeded (i.e. w c < w el because w el is negative), the corresponding part of the foundation enters path 2 . Furthermore, the separation between the rail and the supporting structure is allowed in this constitutive model if the displacement of the beam w c is greater than the plastic deformation w pl (which has a negative value) provided that the plastic deformation has already developed (i.e. w pl = 0); at the locations where no plastic deformation develops, there is permanent contact between the beam and the supporting structure. Besides the foundation stiffness, also the foundation damping is modified if the separation occurs according to the following equation: It must be emphasised that the specific constitutive relation is not of crucial importance for the goal of this paper; other nonlinear constitutive relations can be implemented (e.g. hyperbolic soil model) provided that they are piecewise-linearly defined. The parameters of the constitutive relation are also spatially varying due to the system's inhomogeneity. In reality, the stiffness change is abrupt in some cases and smooth in others. In this paper, the spatial profile of the foundation stiffness is based on a sine squared. The smoothness can be adjusted by changing the transition length l t ; for very small lengths, this profile is close to the piecewise one, for medium lengths it is almost linear, and for large lengths it is smooth. The spatial profile of the foundation stiffness reads [7] where h = {A, B, C, D} denotes the specific branch (Fig. 2); x tc represents the transition centre, and p the stiffness ratio between the stiff and soft domains. To impose a similar spatial variation for the foundation damping, this is defined in terms of the foundation stiffness as is done for a single degree-of-freedom system: where ζ is some damping ratio. Note that, unlike for a single degree-of-freedom system, ζ is not the ratio of actual damping and critical damping. Equations (1)-(10), together with the constitutive relation and the spatial profile of the foundation, constitute a complete description of the problem. In the next section, the solution methods used to solve this problem are described.

Solution methods
As discussed in Introduction, the chosen solution methods are the sequential Laplace transform method, the time-domain method, and the pseudo-force method. To keep the following derivations concise, a bilinear elastic constitutive relation is chosen for the foundation (i.e. only branches A and B from Fig. 2 are considered with the unloading path being the same as the loading one). The extension to more than two branches does not pose any difficulties, and thus, the bilinear law suffices for the demonstration of the following procedures. The bilinear force exerted by the foundation reads where is the stiffness difference between the two branches. Note that w el has a negative value.

Sequential Laplace transform method
The first solution method presented is the so-called sequential Laplace transform method [7]. This method is only applicable when the system's behaviour is piecewise linear, implying that the system behaves linearly between the moments (i.e. nonlinear events) at which its parameters, being functions of the field variables (displacements, velocities, etc.), change abruptly. This enables the application of the Laplace transform over the time between nonlinear events. The solution method described in this section is based on [7], where the approach is described in detail, and thus, only its main aspects are described in the following.
Assuming that the whole system is in the first branch (with stiffness k A d,c ) of the constitutive relation at t = 0 (assumption made in Sect. 2.1), the forward Laplace transform is applied over time to the governing equations, Eqs. (1)- (6). The Laplace-domain solution in the computational domain cannot be obtained analytically for all stiffness and damping profiles. Hence, the computational domain is discretised by means of the finite difference method. The discretised Laplacedomain equation of motion valid in the computational domain reads whereŵ c,1, j is the Laplace-domain displacement vector, I i j is the identity matrix,f IC c,1,i is the forcing vector induced by the initial conditions,f ML c,1,i is the Laplace-domain force exerted by the moving load, and s = σ + iω is the Laplace variable, where σ is a small and positive real number and ω represents the angular frequency. Furthermore,K i j is the bending-stiffness matrix of the beam that includes the contribution of the boundary conditions that is dependent on the unknown displacement, whilef B c,1,i represents the contribution of the boundary conditions independent of the unknown displacement. (This part originates from nonzero initial conditions of the left domain and from the moving load acting on the adjacent domains and is discussed later on in this section where the non-reflective boundary conditions are introduced.) The non-reflective boundary conditions to be applied to the computational domain are derived later in this section; for now, the boundary conditions of the computational domain are kept general. Finally, subscript 1 represents the first time interval. The solutionŵ c,1, j of the linear system of algebraic equations [i.e. Eq. (15)] can be straightforwardly obtained by using a solver for linear systems. To obtain the time-domain solution, the inverse Laplace transform ofŵ c,1, j is numerically evaluated. The resulting time-domain solution is applicable until the first nonlinear event τ 1 (i.e. when w el is first exceeded). The response of the system after the first nonlinear event is governed by an altered equation of motion, which is obtained by changing Eq. (1) (for x ∈ [0, L]) as follows: 1. The reaction of the Kelvin foundation is modified by assigning the adequate stiffness to the nodes where the elastic limit has been exceeded. 2. A time-variable change is introduced, which reads t 2 = t − τ 1 for t ≥ τ 1 (t represents the global time). A graphical representation of the different time axes and nonlinear events is given in Fig. 3. 3. Continuity is ensured by imposing the displacement and velocity of the previous system at τ 1 as initial conditions for the new system (t 1 = t): 4. The boundary conditions (discussed later in this section) as well as the position of the moving load are 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 vari-able t 2 and the procedure described above is repeated. To this end, the procedure is generalised. The Laplacedomain equation of motion for the nth time interval reads where s n represents the Laplace variable associated with the time variable t n . The force exerted by the foundation is split into its contribution (k d,c,n,iŵc,n,i ) proportional to the unknown displacement, and the contribution independent of the unknown displacement, which is accounted for through the external forcê f NL c,n,i , with superscript NL standing for nonlinear. Their expressions are given as follows: Repeating the described procedure until reaching the final time moment t max of the simulation, the displacement vector of the computational domain becomes where t is the time spacing and N is the index of the last time interval.
One important aspect that requires mentioning is that the continuity of displacements and velocities at the nonlinear events is of crucial importance for obtaining correct results, and this depends on the accuracy of the numerical evaluation of the inverse Laplace transform for t n = 0. 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. A method of incorporating the high frequencies without increasing the computational effort is presented in Ref. [7] and is based on asymptotic approximations of the high-frequency behaviour of the two quantities (displacement and velocity).
Up to this point in the derivations, the boundary conditions of the computational domain were kept general.
The derivation of the non-reflective boundary conditions is thoroughly described in Ref. [7]. The forward Laplace transform is applied over time t n to the equations of motion of the two semi-infinite domains, Eq.
(1) for x ∈ (− ∞, 0] and x ∈ [L , ∞), as well as to the interface and boundary conditions, Eqs. (2)-(6). By solving the resulting two boundary-value problems, 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. These reaction forces are the boundary conditions to be imposed on the computational domain and read where the entries in the matrices represent the dynamic stiffness coefficients giving rise to the boundary forces dependent on 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 are given in Eqs. (46) and (47) in Ref. [7]. In addition, vectorsb IC l,n andb IC r,n contain the influence of the initial conditions prescribed to the left and right domains on the reaction forces, giving rise to boundary forces independent of the unknown displacement and slope of the computational domain. Similarly, vectorsb ML l,n andb ML r,n contain the contribution of the moving load to the reaction forces. (They are nonzero only when the moving load is present in the corresponding domain.) The expressions of these vectors for the scenario when x e = 0 are given in Ref. [7], and can be extended for the situations when x e < 0. The expressions are not presented here for brevity.
The non-reflective boundary conditions [Eqs. (21) and (22)] for the computational domain are now fully determined. The contribution of the boundary conditions which is dependent on the yet unknown displacement and slope is incorporated into the beam's bending stiffness matrixK i j , while the contribution which is independent of the unknown displacement and slope is accounted for through the boundary-forcing vector f B c,n,i [see Eq. (18)]. As can be seen, the beam's bend-ing stiffness matrix does not change from one system to the other; however, the boundary-forcing vector needs to be updated at each system change.

Time-domain method
The time-domain method is one of the more conventional approaches to solve the current problem. The spatial dimension is discretised using the finite element method and the Newmark-β time-stepping method is used to solve the discretised system. The only difficulty arises when implementing the non-reflective boundary conditions. As done in the sequential Laplace transform method, the semi-infinite domains are treated analytically. In the time domain, this is done through convolution integrals, which implies that the relations between displacements and forces at the boundaries are history dependent. Therefore, after spatial discretisation, a system of integro-differential equations has to be solved. The approach followed here differs from those in other works using other absorbing boundaries (e.g. [49,50]) in the following ways: (i) in principal, it introduces no numerical reflections (instead of mitigating reflections); (ii) it does not increase the number of degrees of freedom (like in the case of absorbing layers, that model buffer domains with the intention of attenuating waves); (iii) it allows the approach and exit of the load from the modelled domain, thus avoiding transients due to sudden entrance/exit of the load, and in this way it limits the model to the vicinity of the region with support variations and thus further decrease the number of degrees of freedom. The procedure is described in more detail in the following. The equation of motion for the computational domain after the spatial discretisation reads  (14)], respectively. The moving-load forcing vector for a single element f ML c,i is obtained as follows: where φ is the shape-functions vector. The assembly of the global moving-load vector is done in the traditional way and results in a time-dependent vector which has nonzero entries only at the nodes related to the element where the moving load is acting. Similarly, by assuming a constant nonlinearity force inside one element, the nonlinear-forcing vector corresponding to the bilinear constitutive relation is obtained as follows: As for the non-reflective boundary conditions, the expressions derived in Sect. 2.3 are valid for this method too. However, they need to be expressed in the time domain. Moreover, the dynamic stiffness coefficients (k h,Vυ ,k h,Vφ ,k h,Mφ , andk h,Mυ in Eqs. (21) and (22), where h = {l, r}) increase with increasing frequency, and therefore, it is difficult to obtain their counterparts in the time domain. Consequently, for the timedomain method, instead of making use of the dynamic stiffness coefficients (i.e. imposing forces as boundary conditions), the dynamic compliance coefficients are used (i.e. displacement and slope are imposed as boundary conditions), which decay with increasing frequency. By expressing the displacement and slope in terms of the bending moment and shear force in Eqs. (21) and (22), the dynamic compliance coefficients in the Laplace domain are obtained, and their expressions read where and the branches of the complex-valued wavenumbers are chosen such that Im(k h ) < 0 and Re(k h ) > 0. The time-domain non-reflective boundary conditions are now expressed through convolution integrals. Their expressions are given as follows: where w V and w M represent the displacements caused by an applied shear force and displacement caused by an applied bending moment, respectively, while w V and w M represent the same quantities but for the slope. Furthermore, their expressions read . The integrals in Eqs. (28)-(31) are discretised, and the forces are assumed to be constant during one time step and equal to the average between the previous and the next time step [51]. The remainder of the procedure is presented only for w V (0, t) since the procedure is exactly the same for the other terms. The expression for w V now reads where H represents the response of the semi-infinite system to a unit step force applied at t 0 = 0. The response H can be obtained with a good accuracy using the inverse Laplace transform as follows: To obtain accurate results without integrating up to very high frequencies, the asymptotic approximation approach described in [7] is used. Substituting Eqs. (33) and (34) in Eq. (32), and rearranging the terms, the displacement w V becomes It is important to note that in Eq. (35) and the similar expressions for the other three components w M , w V , and w M , the displacement, slope, bending moment, and shear force at time moment t n are unknown, while all the other components are known (i.e. history terms). Therefore, Eq. (35) is divided into a yet unknown instantaneous contribution and an already known history contribution: After deriving the expressions for the other three components (i.e. w M , w V , and w M ) and substituting them in Eq. (27), the non-reflective boundary conditions become: The force and moment are obtained from Eqs. (38) and (39) where vectors b hist l and b hist r incorporate the history shear forces and bending moments and read b hist The forces expressed in Eqs. (40) and (41)  to the Euler-Bernoulli stiffness matrix K EB because they are dependent on the unknown displacement and slope. The remaining terms in Eqs. (40) and (41) are accounted for in the boundary-forcing vector f B c because they are independent of the unknown displacement and slope, and are thus treated as external forces. Incorporating these boundary conditions [Eqs. (40) and (41)] ensures that the finite computational domain behaves as an infinite one. By discretising the convolution integrals, the system of integro-differential equations is approximated by a system of coupled ordinary differential equations with state-dependent coefficients. This system is solved by means of the Newmark-β method [27]. To this end, the generalised displacement vector x n+1 at time moment t n+1 is expressed as a function of the displacement x n , velocityẋ n , and accelerationẍ n at time moment t n as follows: where γ and β are two parameters that indicate how much of the acceleration at the end of a time interval influences the relations for the velocity and displacement at the end of that interval [27]. The parameters are chosen as γ = 1 2 and β = 1 4 , implying that the acceleration is constant over a time step and is equal to the average between the previous and the next time step. This choice is preferred in order to have consistency between the assumed force evolution of the non-reflective boundary conditions [as expressed in Eq. (32)] and the assumed force evolution in the timestepping scheme. After obtaining x n+1 , the generalised velocity and acceleration vectors are also computed to be used for obtaining the generalised displacement at the next time moment. Their expressions reaḋ At each time step, the nodes are monitored to check in which branch of the constitutive relation they are, and if they have changed branch, the Kelvin stiffness matrix K as well as the nonlinear-forcing vector f NL c are updated accordingly. It must be noted that when the nonlinear constitutive relation discussed in the problem statement is adopted, also the Kelvin damping matrix C becomes state-variable dependent and thus needs to be updated at each nonlinear event. The time-domain solution method is now fully described. Next, the pseudoforce method is introduced.

Pseudo-force method
In this section, the problem is solved using the framework of the pseudo-force method [44,45]. In accordance with this approach, the solution method is based on the response of the linear and piecewisehomogeneous system, which in this work is expressed in terms of the Green's function. The steps of the procedure are as follows. The base system is assumed to be linear and piecewise-homogeneous while the nonlinear and the remainder of the inhomogeneity components [the difference between the piecewise-homogeneous profile and the one described by Eq. (12)] of the system are accounted for by means of pseudo-forces. Basically, the nonlinear and the remainder of the inhomogeneity terms are moved to the right-hand side of the equation of motion and the resulting implicit equation is solved in an iterative manner. It must be noted that the shiftinvariant homogeneous system could be chosen as the base system too; however, by doing so, the pseudoforces that account for the inhomogeneity must act on the entire right semi-infinite domain. For the piecewisehomogeneous system, the pseudo-forces need to be prescribed only in the transition zone and its vicinity and this is the reason why the piecewise-homogeneous base system is preferred. The procedure on which this solution method is based is explained in detail in [45]. Therefore, only the crucial aspects of the approach are elaborated in the following. After moving the terms accounting for the nonlinearity and for the remainder of the inhomogeneity to the right-hand side, the equation of motion reads where w is the displacement of the entire beam (i.e. x ∈ are the piecewise-homogeneous damping and stiffness of the Kelvin foundation, respectively, with x tc being the location of the abrupt transition, and f u and f v are the pseudo-forces; f u is proportional to the displacement while f v is proportional to the velocity, and, in the case of the bilinear constitutive relation, they read It must be noted that the part of the nonlinear forcing that is independent of the unknown displacement (denoted by superscript NL in the other two methods), , is incorporated here in f u because in this method both contributions (dependent on or independent of the unknown displacement) are treated as external forces; therefore, the distinction between them is not needed in this method. Additionally, it must be noted that when the nonlinear constitutive relation discussed in the problem statement (Sect. 2.1) is adopted, the behaviour of the foundation damping becomes piecewise linear too and must be accounted for in Eq. (51). The solution to Eq. (49) can be expressed as a superposition of the response w ML caused by the moving load f ML in the base system and the response to the pseudo-forces. Firstly, the response w ML caused by the moving load can be obtained 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 (e.g. [23,46]). To obtain the time-domain solution, the inverse Fourier integral is evaluated numerically. Secondly, the response to the pseudo-forces is expressed as a convolution integral of the impulse response of the base system and said forces. For conciseness, the derivation is demonstrated for f u only. The contribution w v caused by the forcing term related to the damping inhomogeneity f v is kept general and is made specific in the final expression [Eq.
(64)]. The response thus reads Fig. 4 The assumed spatial profile of the pseudo-forces; the length of the domain should be chosen such that the first (i = 1) and last (i = N x ) pseudo-force are nearly zero where g is the time-domain Green's function of the base system. To evaluate Eq. (52), the integrals are discretised. Firstly, the displacement w is assumed to be piecewise constant in space and equal to the value at the centre of each discrete element. Consequently, also the Green's function g is determined with a box-car shaped load in space (Fig. 4), while it assumes a Dirac delta function loading in time.
The continuous-in-time and discrete-in-space expression for the displacement reads where i is the index for the observation point in space whileī is the index for the running (integration) spatial variable, and N x is the number of integration points. Secondly, the forces f ū i are assumed to be piecewise linear in time. The displacement thus becomes where n is the index for the observation time variable whilen is the index for the running (integration) time  54) is valid only for n ≥ 1 because at t 0 = 0, just the response caused by the moving load (i.e. w i,0 = w ML i,0 ) is present, similar to the initial conditions imposed in the other two methods. It can be observed that Eq. (54) consists of two terms, one proportional to the force at time moment tn −1 and one proportional to the force at time moment tn. Since the two forcing terms are not dependent on the variable of integration τ , they can be taken out of the integral. The equation can therefore be rewritten as follows: where L and R represent the responses observed at t n due to triangular pulses lasting between tn −1 and tn (Fig. 5). The time-domain Green's functions g i,ī (t n − τ ) need to be obtained numerically while the integration from tn −1 to tn needs to be performed numerically too, thus introducing two sources of inaccuracy. To increase accuracy, the response associated with a triangular pulse can be obtained directly from the Laplace domain, where only the inverse Laplace transform needs to be evaluated numerically. To this end, the expressions for L and R can be rewritten by introducing the variable changeτ = τ − tn −1 : where t n−n+1 = t n − tn −1 . These responses can now be expressed as inverse Laplace transforms as follows: whereĝ i,ī is the Laplace-domain Green's function associated with a Dirac delta load in time and a box-car function in space, andP L andP R are given bŷ The inverse Laplace transforms are evaluated numerically using a quadratic, nested, adaptive integration scheme.
To express the contribution from the forces proportional to the velocity f v [Eq. (51)], one needs to derive an expression for the velocityẇ i,n . One could follow the same procedure as for the displacement w i,n and obtain a similar equation to Eq. (55). Then, the equations for w i,n [Eq. (55)] andẇ i,n could be solved simultaneously. However, the Laplace-domain counterparts of L and R that relate applied force to resulting velocity have a poor decay and, consequently, evaluating the inverse Laplace transforms is time consuming. A computationally efficient alternative is to approximate the velocityẇ i,n as a function of the displacement at the previous time moments by using the finite difference method. This makes the force f v i,n proportional to the displacement (i.e. f v i,n (w i,n , w i,n−1 , w i,n−2 )), so that it can be added to the force f u , together being incorporated in the force f u+v i,n . Note that for conciseness, the force f u+v i,n is indicated in the following to be dependent on the displacement w i,n , while in actual fact it is also dependent on the displacement at the previous time moments. The complete solution now reads It can be observed that Eq. (64) is implicit forn = n. Therefore, the equation is divided in an yet unknown instantaneous contribution and an already known history contribution. The equation becomes In order to advance to the next time step, Eq. (65) is solved for w i,n using an iterative scheme. The scheme is defined by the following recursive relation that starts at j = 0 (where j indicates the iteration counter): The rate of convergence of this iterative scheme is R i,ī ,0 [45]. Iterations continue until specified tolerances are met for all entries of the displacement w i,n . No convergence problems were encountered while computing the results for this paper. All three solution methods have been presented now. In the following section, the performance of the three solution methods is compared.

Results and discussion
Here, the proposed solution methods, namely the sequential Laplace transform (SLT) method, the timedomain (TD) method, and the pseudo-force (PF) method, are compared in terms of accuracy, computational efficiency, and feasibility to apply the methods to more complex models. The three methods are assessed for extreme situations. To this end, a relatively high load velocity is chosen, namely 95% of the minimum phase velocity in the soft domain (the velocity is 400 m/s, which may seem extremely high; however, the velocity relative to the minimum phase velocity is of importance, and not the absolute value of the velocity; for a more extensive discussion on this subject, see Sect. 3.1. of Ref. [26]). Also, a relatively low damping ratio (ζ = 0.05) is used to make sure that errors are not damped out. Moreover, a high stiffness dissimilarity between the left and right domain ( p = 5) is chosen which, combined with an abrupt transition (l t = 0.01 m), leads to strong transition radiation and large plastic deformation. The location of the abrupt transition is at x tc = 5 m for all simulations presented in this section. Also, x e = − 30 m for all the results presented in the following to ensure that the eigenfield does not interact with the inhomogeneity at the beginning of the simulation. The parameter values used in the computations are given in Table 1.

Accuracy
The accuracy assessment of the three methods is based on two types of simulations. Firstly, a linear limit case is considered by setting the parameters k B d,c (x), k C d,c (x), and k D d,c (x) to be equal to the initial stiffness k A d,c (x). By doing so, fictitious nonlinear events are introduced, but the response should be the same as the response of the linear system. In the second type of simulation, the nonlinear system with the constitutive relation from Sect. 2.1 (Fig. 2) is considered.
The linear limit case is compared to the semi-analytical solution of a piecewise-homogeneous and linear system. This semi-analytical benchmark solution is the same as the response w ML of the base system to the moving load in the PF method (Sect. 2.5). Because the only source of error of this semi-analytical solution is the numerical evaluation of the inverse Fourier transform, a very high maximum frequency (i.e. f max = 40 kHz) and a small frequency step (i.e. f = 0.1 Hz) are adopted to obtain a very high accuracy. The relative error e(x) used in this section is defined as follows: where w lin is the displacement in the linear limit case and w bench is the benchmark solution. This error formulation where the summation over time is performed in both the numerator and denominator is chosen because at certain locations and time moments, the displacement can be zero or close to zero, which, if the summation over time were not taken, would lead to a huge relative error that has no physical significance. This error relation is also used for the nonlinear case; however, for the nonlinear simulation there is no semi-analytical solution to use as a benchmark, and thus, the three methods are compared to each other. The accuracy is studied for varying sampling parameters t and x. The three methods have different sensitivities to the two sampling parameters; therefore, when varying one sampling parameter, the other one is chosen such that the resulting error is of the same order for all three methods. For the SLT and TD methods, the frequency spacing ω = 4π rad/s has been chosen after sensitivity studies, while for the PF method the adaptive algorithm (see Ref. [45]) chooses itself the frequency sampling. Note that for the TD method, the frequency sampling is needed for evaluating the inverse Laplace transforms to obtain the response functions H [see Eq. (34)] and the boundary vectors (b IC l , b ML l , and b ML r ). As for the maximum frequency, f max = 2 kHz for the TD method, f max = 1/(2 t) for the SLT method, and f max = 2 log 2 (2π/ t)+2 for the PF method, unless specified otherwise. It must be mentioned that f max for the PF method is set to such a large value (e.g. f max ≈ 40 kHz for t = 100 μs) because the adaptive algorithm chooses very few sampling points if the integrand is smooth at high frequencies, therefore leading to almost no additional computational effort.

Linear limit case
The elastic displacement limit ratio q = w el /w e max is chosen to have a value close to unity (q = 1.01) such that fictitious nonlinear events occur in abundance. For the TD method, this limit case is the same as a linear case (no fictitious nonlinear events) in terms of accuracy. For the SLT method, this limit case tests all operations, namely solving the system in Eq. (18), updating the non-reflective boundaries, and evaluating the inverse Laplace transform. The PF method reduces to computing the response w ML caused by the moving load, which is exactly the same as the benchmark solution; therefore, such a comparison would be superfluous. To also test the PF method, the piecewisehomogeneous base system is prescribed to have the abrupt transition at x tc = 6 m, while the forces that account for the inhomogeneity need to act such that the system simulates an abrupt transition at x tc = 5 m.  In this way, the Green's functions and the convolutions are assessed. Figure 6 shows that the response obtained with all three methods converges to the semi-analytical one as x decreases. However, to obtain similar magnitudes of the error, the SLT method requires a smaller x than the TD and PF methods. This is because a certain ratio between the maximum frequency f max and x must be satisfied in the SLT method. More specifically, the Laplace-domain moving load is harmonic in space ; to accurately representF ML at high frequencies, there should be at least five spatial points per wavelength, leading to the requirement that x = v/(4 f max ). For the parameters in Fig. 6, this minimum requirement is x = 5 cm; for this reason, the error is much larger for x = 10 cm while for the remaining values of x the error is of similar magnitude. The larger error in the right domain for the SLT and TD methods is caused by the significantly smaller values of the displacement encountered there. The PF method is not so much affected by this because as the transients caused by the transition die out, the solution tends to the semi-analytical baseline response.
The error also decreases with decreasing t for all three methods, depicted in Fig. 7. Unlike for decreasing x, here the TD and PF methods require smaller t than the SLT method to obtain a similar error magnitude. The TD method exhibits a significant error reduction with decreasing t, which is to be expected since the Newmark method is a time-stepping method and the smaller the t, the more accurate the approximation becomes. For the SLT method, f max needs to be increased with decreasing t to obtain convergent results. More specifically, to satisfy the continuity conditions [see Eqs. (16) and (17)], the inverse Laplace transform requires a high accuracy close to t n = 0 (here, n refers to the time interval; see Sect. 2.3). This can be achieved by satisfying the criterion t = 1/(2 f max ) (i.e. the time discretisation should not be more refined than the Nyquist criterion requires). The fact that the spatial discretisation is linked to f max (as stated in the previous paragraph) leads to the three sampling parameters (i.e. x, t, and f max ) being interdependent for the SLT method. As presented in Sect. 3.2, this leads to significant increase in computational effort when one of them is changed. As for the PF method, the error decreases considerably with decreasing t. This is caused by the pseudo-forces (acting between x = 5-6 m) becoming more accurate; in other words, the approximation of the convolution integrals [see Eq.

Nonlinear case
As there is no semi-analytical solution to test against in the nonlinear case, it is chosen to run a high-resolution simulation using the TD method (i.e. f max = 4 kHz, x = 0.5 cm, t = 10 μs), and perform the accuracy analysis using this result as a benchmark. Therefore, it is important to realise that the error studies presented in this subsection do not represent errors with respect to a true solution, but they are essentially differences between solution methods. For the results presented in this subsection, the elastic displacement limit ratio is chosen as q = 1.05.
One can observe at a first glance of Figs. 8 and 9 that the error in general is considerably larger than for the linear limit case (Sect. 3.1.1). The additional error is mainly caused by two factors: the spatial and temporal discretisation of the foundation's nonlinear force [see Eq. (14) for the bilinear case], and the combination of non-smooth nonlinearity and separation between the beam and supporting structure only being allowed in the locations where plastic deformation has developed (see Fig. 2). To exemplify the latter, assume the minimum displacement at a certain location x 1 computed with the TD method to be min w(x 1 , t) = w el + ε (where ε is a very small displacement); this location is therefore in the initial stiffness branch (k A d,c ) for all time moments. Now, assume the same quantity computed with the SLT method to be min w(x 1 , t) = w el − ε; this location experiences a small plastic deformation and the separation between the beam and foundation is allowed (i.e. no force is exerted by the foundation for upwards displacement of the beam). As can be seen, a very small difference in the displacement magnitude (2ε) can lead to different behaviour, causing a larger error than in the linear limit case. Nonetheless, the plastic deformation profile and the displacement's time history show that despite the substantial errors, the solutions exhibit the same trend for all three methods (Figs. 10, 11). Figure 8 shows that all solutions converge to the benchmark by decreasing x; however, the convergence is not monotonic for all locations, not even for the TD method despite the benchmark solution being computed with the TD method. The non-monotonic convergence is caused by the several sources of error that may eliminate each other, leading to apparently accurate solutions for certain sampling combinations. Interestingly, the largest error generally does not occur where the plastic deformation develops, but at a loca-tion (around x = 3.5 m) between two plastically deformed regions (Fig. 10).
One can observe in Fig. 9 that the error decreases with decreasing t for all three methods, but the decrease is non-monotonic similarly to the error for decreasing x, and can be explained in the same manner. The errors obtained with the TD and PF methods seem to converge to a non-zero value because with the decrease in t the error introduced by the spatial discretisation starts governing. This is not observed in the SLT method because x is decreased with t to respect the criterion To offer a bigger picture, Figs. 10 and 11 present the plastic deformation profile and the displacement's time history at x = 3.5 m (the location in the soft domain with the largest error), respectively. Although in some cases the errors may appear to be large (e.g. 2-20% in Figs. 8,9), all three methods perform very well and the differences can hardly be seen.
To conclude, all three methods converge given the correct sampling parameters x and t, and all three can reach similar error magnitudes. Nonetheless, also the computational effort required to achieve the observed accuracy needs to be investigated if one wants to present a complete picture. In this section, the sampling parameters have been chosen such that the error magnitude is similar for the three methods. In the next section, the computational effort to obtain the accuracy discussed in this section is presented and discussed.

Computational efficiency
In this subsection, the most intensive computational operations are identified for each method. These computational operations are analysed qualitatively as well as quantitatively by presenting a comparison of the computational time required to perform the simulations presented in Sect. 3.1.2. The computational effort presented in Figs. 10 and 11 represents wall-clock time and the simulations have been performed on a PC with Intel(R) Xeon(R) CPU E5-1630 v3 @ 3.70 GHz processor and 32 GB of RAM memory. It must be emphasised that further optimisation of the algorithms is always possible; however, the authors made their best to optimise them and believe that further improvement, although possible, will not lead to significantly different results, especially for the qualitative analysis.
For the TD method, the most computationally demanding processes are numerically solving the system of Eqs. (44)- (46) and computing the non-reflective boundary conditions. Firstly, determining the solution of the system requires the computation of the matrix inverse Y −1 n+1 (the most computationally intensive operation) once at the start of the simulation and recomputed at each nonlinear event due to the change in the properties of the foundation (i.e. K and C). Alternatively, one can solve the linear system of equations Y n+1 x n+1 = z using a sparse Y matrix; this operation is considerably faster than explicitly computing Y −1 . However, this operation would need to be computed for each time moment while the matrix inversion needs to be computed only at each nonlinear event. Consequently, for few nonlinear events, the former option is faster while for many nonlinear events, the latter is more efficient. Secondly, the computation of the non-reflective boundary conditions involves numerically evaluating the inverse Laplace transforms to obtain the response functions H [see Eq. (34)] and the boundary vectors (b IC l , b ML l , and b ML r ), and evaluating the convolution integrals at the boundaries [see Eq. (32)].
For the SLT method, the most computationally intensive operations are solving the linear system of equations [Eq. (18)], updating of the non-reflective boundary conditions, and numerically evaluating of the inverse Laplace transform. The linear system of equations [Eq. (18)] needs to be solved for each frequency and at each nonlinear event. This is done by using an algorithm to solve linear systems of equations in MAT-LAB (i.e. the mldivide routine) where the dynamic stiffness matrix is assembled as sparse. For the results presented here, this is the most computationally intensive part of the SLT method. Next, updating the boundary conditions involves computing the particular solutions of the non-trivial initial conditions in the left and right domains; this implies a numerical evaluation of the inverse Laplace transform for each frequency (see Eqs. (68) and (69) of [7]). Finally, the numerical evaluation of the inverse Laplace transform needs to be performed for each time moment and for each spatial point.  Figure 12 presents the computational effort C (measured in seconds) needed to reach the accuracy shown in Fig. 8. The TD and SLT methods require similar computational times while the PF method appears to be more computationally intensive. This is caused by the fact that with decreasing x (i.e. increasing the number N x of excitation ξ¯i and observation x i points), the number of required response functions L and R increases with N 2 x . However, once the response functions are determined, a whole parametric study can be conducted (e.g. varying v, l t , types of nonlinearity, etc.) Δt (s) For the TD and SLT methods, the computational effort is governed by the matrix inversion operation while updating the non-reflective boundary conditions, as one could expect, is not affected by the change in x since they need to be determined only at the boundaries. For the SLT method, the inverse Laplace transform appears to be the least computationally intensive operation. Figure 13 presents the computational effort C needed to reach the accuracy shown in Fig. 9. The computational effort in the TD method is still governed by the inversion of the matrix operation despite x being constant because the number of nonlinear events increases with the decreasing t (i.e. the number of matrix-inversion operations increases). Also, the computational effort needed to update the nonreflective boundaries appears to be slightly affected by the decrease in t. Although the evaluation of the convolution integrals increases almost quadratically, the numerical evaluation of the inverse Laplace transforms increases linearly with decreasing t; for t = 80-250 μs, the inverse transforms are governing C NRBC for the TD method, while for smaller t the convolution integrals start governing. As for the SLT method, the matrix inversion is also still governing the computational time when decreasing t because x and f max are also varied for these simulations and the number of nonlinear events increases with decreasing t. Respecting the convergence criteria for the SLT method leads to a considerable increase in computational effort. For the PF method, C appears to increase less drastically than the other two methods with decreasing t. This is because C in the PF method is mainly governed by x, which is kept constant here. Also, C for the PF method is not affected by the number of nonlinear events.
Finally, given the accuracy levels presented in Sect. 3.1.2, the TD method is the most computationally efficient method. However, for certain combinations of t and x, the SLT method is not far behind, while the PF method shows potential for a large number of simulations provided that the associated response functions do not change.
To conclude, it must be mentioned that the abovedrawn conclusions can be somewhat different for another choice of parameters (e.g. train velocity, transition length, elastic displacement limit ratio q, damping ratio, etc.). It is impossible to cover all combinations of parameters, but an indication of how changing these parameters may influence the computational effort is given in the following. The main effects of changing these parameters on the computational effort comparison can be grouped as follows: changes leading to a variation in the number of spatial nodes, changes leading to a variation in the number of time samples, and changes leading to a variation in the number of nonlinear events. Firstly, the number of spatial nodes is affected by the largest real part of the wavenumber of the eigenfield (see Appendix A of Ref. [7]) and of the generated waves; this is in turn influenced by the load velocity. (The higher the load velocity, the larger the real part of these wavenumbers.) The number of spatial nodes is also affected by the length of the transition zone because the complete nonlinearity and inhomogeneity needs to be incorporated in the computational domain. Secondly, the number of time samples is affected by the maximum frequency of the generated waves (same reasoning as for the maximum wavenumber can be applied here), by the load velocity (i.e. the lower the velocity, the more time moments until the load has passed the transition), and by the damping present in the system (i.e. the higher the damping, the fewer time moments until the waves die out in the computational domain). Thirdly, the number of nonlinear events is affected by the load velocity (the lower the velocity, the less plastic deformation and the fewer nonlinear events), by the length of the transition and the stiffness ratio p (the smoother the transition the fewer nonlinear events), and by the elastic displacement limit ratio q (the larger the ratio the fewer nonlinear events). From the previously drawn conclusions, it can be seen that the PF method is most influenced by the variation in the number of spatial nodes, the TD is most affected by the variation in the number of time samples, and the SLT method is most affected by the number of nonlinear events. (It may seem that also the SLT method is strongly affected by the number of time samples, from Fig. 13, but it is mostly the decrease in t which leads to an increase in f max and a decrease in x that con-tributed to this sensitivity; an increase in number of time samples while keeping the t fixed does not follow a similar trend in computational effort increase.) These qualitative observations can give an indication to which method is advantageous for a different choice of parameters.

Feasibility of application to more complex systems
Here, the feasibility for the methods to handle frequency-dependent properties of the structure (e.g. frequency-dependent springs and dashpots) is qualitatively assessed, as well as to handle a smooth nonlinearity (as opposed to the piecewise-linear one), and to apply the solution methods to multi-dimensional models. It must be emphasised that all three methods can be applied in the above-mentioned situations, but one can be more efficient than others, and this is discussed in the following.
Usually when a model reduction is opted for (e.g. going from 2-D to 1-D), the conversion of parameters is done for a specific frequency range of interest. If more frequency ranges are of interest or if the frequency range is broad, then the conversion may lead to frequency-dependent parameters of the reducedorder model. The TD method is inefficient in this case because additional convolution integrals are introduced that are distributed throughout the spatial domain (opposed to the convolution integrals at the boundaries which are local), leading to a significant increase in computational time. For the SLT method, the frequency-dependent parameters can directly be incorporated in the equation of motion with no additional effort since the system is solved in the Laplace domain; therefore, the computational effort in this case is unchanged. For the PF method, the response functions L and R are determined from the Laplace domain and can easily incorporate frequency-dependent parameters too; however, the fact that the nonlinearity and part of the inhomogeneity is accounted through the pseudo-forces leads to additional convolution integrals throughout the spatial domain, which increases the computational effort of the method. Therefore, the SLT method appears to be most suitable to tackle a system with frequency-dependent properties.
In the problem statement, the reaction of the Kelvin foundation was assumed to be piecewise linear. Although under certain conditions (see [7]) this assumption is rea-sonable, other analyses might require a smoothly nonlinear constitutive law (e.g. cubic super-linear model [52]). The TD method can easily incorporate such behaviour, but with an increase in the computational time because the matrix-inversion operation [see Eq. (44)] must be performed at each time moment; however, this operation can be sped up as discussed in Sect. 3.2. The SLT method is very inefficient for such behaviour; the smoothly nonlinear behaviour implies that nonlinear events occur essentially at each time moment leading to a considerable increase in the computational time. The PF method, however, can easily handle such a behaviour. The only increase in the computational time might come from the poorer convergence of the iterative scheme; however, this is unlikely. To present a more quantitative comparison, fictitious nonlinear events are imposed at each time step for one set of sampling parameters ( x = 5cm and t = 80 μs for the TD method, x = 5 cm and t = 250 μs for the SLT method, and x = 5 cm and t = 100 μs for the PF method). As can be seen in Fig. 12, the computational effort for these sampling parameters in the reference case is approximately 40 s for the TD method, 100 s for the SLT method, and 2,660 s for the PF method. When the fictitious nonlinear events are imposed at each time moment, the computational effort is approximately 55 s for the TD method and 550 s for the SLT method, while the one of the PF method remains approximately the same. Thus, by implementing a smoothly nonlinear constitutive law, the TD appears to be most efficient and the PF least efficient, while the computational effort of the SLT method is most affected by such a change.
As for using the methods to analyse a 2-D model, for example, all three methods have advantages and disadvantages. If the nonlinearity and inhomogeneity of the system is considered throughout the depth (i.e. ballast layers as well as the soil layers), then all three methods need to discretise a considerable portion of the vertical direction, leading to a significant increase in the computational time. For such a system, the TD method is likely to be most efficient as the SLT and PF methods become unfeasible for very large systems. More specifically, for the SLT method, the matrix inversion for systems with a very large number of degrees-offreedom becomes unfeasible when performed for each frequency and each nonlinear event, while for the PF method, the large number of associated response functions to be determined in the time domain renders the method unfeasible. However, if the nonlinearity and inhomogeneity are limited to the surface layers (e.g. nonlinear ballast layer resting on inhomogeneous soil layer supported by a linear half-space), then only the top layers need to be discretised and the others could be included by using frequency-dependent parameters. For such a system, the TD method requires additional spatially distributed convolution integrals to account for the presence of the half-space. The SLT method can exactly replace the half-space by frequency-dependent springs and dashpots, which does not affect its efficiency. Also the PF method can incorporate the halfspace with no additional effort (it will already be incorporated in L and R, which are computed only for the top layers). An example of such a 2-D system has been formulated in Ref. [17] (but with a linear behaviour for all materials), where the ballast and the foundation are modelled as lattices. This model has approximately 60,000 degrees-of-freedom (DOFs) and 40,000 of them are for the foundation. Replacing the foundation by frequency-dependent springs can reduce significantly the amount of DOFs; however, for the SLT method, the dynamic stiffness matrix for the reduced system needs to be inverted for each frequency individually, while for the PF method, obtaining the Green's functions even for this reduced system seems unfeasible. To conclude, even though the SLT appears to have the potential to be more efficient than the TD method for such a system, it difficult to draw a definite conclusion.

Conclusions
This study presents three solution methods to obtain the response of a 1-D model consisting of an infinite Euler-Bernoulli beam resting on a locally inhomogeneous and nonlinear supporting structure subject to a moving constant load. The three methods, namely the sequential Laplace transform (SLT) method, the time-domain (TD) method, and the pseudo-force (PF) method, were chosen from the three main categories of solution methods available in the literature, namely integraltransform methods, time-domain methods, and hybrid methods, respectively. The three methods were compared in terms of accuracy, computational efficiency, and feasibility of application to more complex systems.
Results show that all three methods are able to reach similar accuracy levels, both for a linear limit case and for the nonlinear situation. For the latter case, the combination of non-smooth nonlinearity and separation between the beam and supporting structure only being allowed at the locations where plastic deformation has developed leads to larger relative errors than in the former case and to non-monotonic convergence with decreasing x and t. Given similar accuracy levels, the TD method is overall the most computationally efficient method. However, for certain sampling combinations, the SLT method is not far behind while the PF method shows potential for a large number of simulations provided that the response functions do not change. As for the feasibility to apply these methods to more complex systems, the SLT method appears to be most efficient in dealing with frequency-dependent parameters while the TD and PF methods in dealing with smooth nonlinearity. For a 2-D system, if the nonlinearity and inhomogeneity is considered throughout the depth, the TD method is likely to be most efficient; however, if the nonlinearity and inhomogeneity are limited to the surface layers (e.g. nonlinear ballast layer resting on inhomogeneous soil layer supported by a linear half-space), the SLT method has the potential to be more efficient than the TD method.
Based on this study, the adequate solution method to solve more complex systems can be selected depending on the specific characteristics of the model. Although the 1-D model presented in this study has mainly been used to assess the three methods, it can also be used for preliminary designs of transition zones in railway tracks. The part of the piecewiselinear foundation force that is independent of the unknown displacement k l,Vυ ,k l,Vφ ,k l,Mυ ,k l,Mφ Entries in the dynamic stiffness matrix of the left semiinfinite domain k r,Vυ ,k r,Vφ ,k r,Mυ ,k rs,Mφ Entries in the dynamic stiffness matrix of the right semi-infinite domain K i j Bending-stiffness matrix of the beam τ n

Nomenclature
The moment of the nth nonlinear event in the global time axis t τ n The moment of the nth nonlinear event in the local time axis t n I i j Identity matrix N The index of the last time interval s n The Laplace variable associated with the time interval t n t n The nth time interval (Fig. 3) TD method