A lattice model for transition zones in ballasted railway tracks

The procedure presented in this work aims at providing a framework for studying set- tlement of ballast in zones with stiffness variation of the railway track support. The proposed procedure results from expanding an existing inﬁnite periodic model of a railway track to account for variations in the stiffness of the foundation. Ballast is simulated via a linear lattice, whose dynamic response differs from that of a continuum. The expanded model is composed of three regions: a left region, which is semi-inﬁnite and periodic; a mid-region, of ﬁnite length and where the properties of the foundation can change; and a right region, which is also semi-inﬁnite and periodic and whose properties can differ from those of the left region. The equations of motion of the mid-region are written directly in the time domain, with the rail being described by ﬁnite elements. On the other hand, the left and right semi-inﬁnite regions are treated semi-analytically in the frequency domain, and afterwards their responses are converted to the time domain, resulting in convolu-tion integrals prescribed at the boundaries of the mid-region that simulate non-reﬂective boundaries. The ﬁnal model only contains the degrees of freedom corresponding to the mid-region (which can be as short as the region where the stiffness variation is present), and that leads to faster calculations than if the boundaries were placed further away to dissipate undesired reﬂections. The method is cast in the time domain, and all elements are assumed to behave linearly. In the future, the model will be expanded to incorporate non-linear behaviour of the ballast. The presented method is validated by means of simple examples, and afterwards applied to a real scenario in which a culvert crosses a railway track. As presented, the method can be used to study the linear dynamics of transitions zones, study mitigation measures, and infer about indicators like force transmitted and energy dissipated, which might be useful to assess the development of settlement of the ballast.


Introduction
Previous studies and measurements indicate that differential ballast settlement occurs at a faster pace in regions where railway tracks experience variations in their support stiffness than in regions with homogeneous support [1 , 2] . The accelerated differential settlement demands for more frequent maintenance operations, which translates into additional costs and disruptions of the train schedules. These regions of stiffness variation are therefore of interest to railway researchers, since understanding the phenomena behind the accelerated degradation helps designing solutions to mitigate the issue.
Studies focusing on settlement of ballast can be categorized into field measurements and condition monitoring [ 3 , 4 ], in which the position of rails and the geometry of ballast are measured and monitored in the period between maintenance operations; lab tests [5][6][7] , in which ballast layers are cyclically loaded and the evolution of its geometry is recorded; and numerical simulations, in which more or less complex mathematical models are formulated to predict the response of the track. Whilst the first category is very useful for detecting main features of the settlement, it is less viable to optimize mitigation measures, since any alteration would require new constructions, and that would first require an approval from the relevant safety authorities, which are usually very conservative and strict. On the other hand, lab tests allow for a better control of the loading process and in this way are useful to derive constitutive laws for the ballast (which can include mitigation measures such as geogrids [8] or binding polymers [9] ), but for studying railway track support stiffness changes this approach becomes unfeasible due to the space required to accommodate such features. Numerical simulations have the disadvantage of belonging to a fantasy world, in which everything is well defined and controlled, but are very versatile and thus useful for assessing the efficiency of mitigation measures. The three categories of studies complement each other. In this work, the focus lies on numerical simulations.
A wide variety of models focussing on longitudinal variations of the track support have been presented in the literature. They include simple 1D models of rails resting on viscoelastic Winkler foundations with longitudinally varying properties [10][11][12][13] , simple 2D models of half-spaces with step variations in their properties [ 14 , 15 ], and very detailed 2D or 3D models of track, foundation and potential mitigation measures [16][17][18][19] . For the latter, the most common solution strategy is the use of the finite element method (FEM) to discretize the domain and its different components, which eventually can be described with non-linear constitutive laws that define the settlement behaviour of ballast [20] . However, because FEM can only handle finite domains, in order to avoid boundary effects, lar ge volumes of domain have to be modelled, leading to huge computational times and rendering the approach not suitable for large parametric studies.
In what concerns ballast description, one can assume continuum models or discrete models. The first assumption is readily available in any FEM software, and can be combined with non-linear constitutive laws that define the ballast behaviour [21] . However, representing ballast in such a way fails to capture certain features of its dynamic response, namely the existence of a threshold frequency above which waves cannot propagate, as well as the specific dispersive behaviour at high frequencies; this may very well affect the local dynamic behaviour of the ballast and thus the development of permanent deformations. Discrete models, on the other hand, intend to describe each particle of ballast independently, and thus can provide more accurate descriptions of its behaviour. The discrete element method (DEM) [22][23][24][25][26] is the method used the most, since it can describe the contact between particles quite well and capture different failure modes of ballast, but due to its high computational cost, it cannot be used to simulate the lengths of the track needed to consider changes in the stiffness of the support. Alternatively to DEM, lattice models [ 27 , 28 ] also intend to describe each ballast particle individually, but they assume that all particles are equal and regularly distributed and that contacts between particles do not change, in this way avoiding the contact calculations and reducing the computational time. By making the lattice connections non-linear, settlement behaviour can be simulated.
The procedure presented in this work aims at providing a framework for studying train-passage induced settlement of ballast in zones with stiffness variations of the support. The purpose of the framework is to allow parametric studies on mitigation solutions, and therefore the calculations should not be too computationally intensive. Therefore, a time-domain 2D model is proposed, in which ballast and foundation are described by a linear lattice (in future work, non-linear lattice connections will be added to look at settlement behaviour), and the infinite character of the track in the longitudinal direction is achieved with non-reflective boundaries. The proposed approach consists in dividing the ballasted track in three regions: one before the stiffness variation that is semi-infinite in the longitudinal direction; one after the stiffness variation that is also semi-infinite; and a middle region of finite length where the stiffness variation is described. Because the domains before and after the stiffness variation are replaced by transmitting boundaries, the final model only contains the degrees of freedom corresponding to the middle region (which can be as short as the region where the stiffness variation is present), and that leads to faster calculations than if the boundaries were placed further away to dissipate undesired reflections. The resulting approach is an expansion of the model presented in a previous work [29] , in which a semi-analytical solution for the track without stiffness variation is presented.
This manuscript is organized as follows. Section 2 introduces the model in terms of the three regions mentioned in the previous paragraph and formulates the equations of motion of each region in the time domain. In Section 3 , the equations derived in a previous work [29] are worked out in order to obtain the transmitting boundaries that replace the semi-infinite regions. In Section 4 , the method is verified by comparing each intermediate step with equivalent models, and the threeregion model with the analytical expressions from [29] . In Section 5 , the usage of the proposed approach is illustrated based on a real stretch of the Dutch railway line. In the last section, Section 6 , the work is summarized and final considerations are proffered.

Three region model
The procedure presented in this article builds on the model introduced in a previous work [29] , in which the railway track is considered to be periodic and is solved semi-analytically for moving and non-moving loads. There, rails are modelled   as an Euler-Bernoulli beam, sleepers as rigid masses, and ballast and foundation as a horizontally layered lattice, which is a regular network of equally spaced masses connected via viscoelastic connections, as defined in the work by Suiker and collaborators [ 27 , 28 ]. The weight of the lattice masses and the properties of the viscoelastic connections may change from layer to layer in order to distinguish the different material properties of ballast and foundation. For the ballast layer, the properties can be inferred from lab compression tests (see Section 4 ), while for soil/foundation layers, the properties can be inferred from equivalent continuous models [27] (be aware that lattices limit the approximation of continuum layers to Poisson's ratio below 0.25, see Section 5.1 ). Fig. 1 schematizes the periodic model.
The expansion described in this work consists in adding a mid-region in which the properties of the foundation (or of any other component, though here the interest is in the foundation) can vary. In this way, the model is divided into three parts: a left region, which is periodic and semi-infinite and formulated semi-analytically in the frequency domain; a midregion, within which material properties can change and that is formulated in the time domain; a right region, which like the left region is semi-infinite and formulated semi-analytically in the frequency domain. The left and right regions may have distinct ballast or foundation properties. Fig. 2 schematizes the division of the sub-regions.
In the following subsections each of the sub-regions is described, and the procedure to couple them is explained. The mid-region is described first. Fig. 3 shows the mid-region isolated from the rest of the domain. In the figure, the mid-region is represented by 3 regular cells of the same type as those of the left and right regions, but in principle the mid-region can consist of any type of structure, provided that the geometries at its left and right boundaries are compatible with the side regions. In addition, the material properties of the mid-region can change from location to location, and its length can assume whatever value necessary. This region is acted upon by external forces ( f ext , not represented in Fig. 3 ), which in most cases will be a set of loads moving on top of the rail from one edge to the other, and by the interaction forces at the left ( f l ) and right ( f r ) boundaries. These interaction forces represent the effect of the semi-infinite regions on the middle domain, and will ensure that the simulation is not affected by unwanted reflections at these edges. The interaction forces on each side are composed by a shear force and a moment at the rail, and a couple of horizontal and vertical forces at each mass of the lattice. Due to its simplicity, the Finite Element method is used to describe the behaviour of the rail, which in this work is divided into elements that are half the length of the lattice characteristic distance d (in principle, any discretization length can be used, but for convenience associated with meshing, d/ 2 is used). Thus, after discretizing the rail and assembling all rail elements and viscoelastic connections of the lattice into the stiffness and damping matrices ( K and C ), and all the inertial elements into the mass matrix M , the equations of motion of the mid-region can be written as

Mid-region
For convenience, the degrees of freedom have been divided into those at the left interface ( l ), those at the right interface ( r ), and the remaining ones ( m ). The sub-indices in the matrices K , C and M , and in time dependant vectors u and f represent these groups of degrees of freedom. In Eq. (1) , the time dependant vector u contains the response of each node (for nodes of the rails and sleepers, the responses are the vertical translation and the rotation; for the lattice nodes, they are the horizontal and the vertical translations). The mass matrices M ll and M rr are zero everywhere except at entries corresponding to the degrees of freedom associated with the rail (note that at the lateral boundaries of the mid-region there are no lattice particles; these are considered in the left and right regions; see Fig. 4 for the definition of the left and right regions). Also, the mass matrices M lm = M T ml and M rm = M T mr are non-zero because the application of the finite element method to describe the rail leads to inertial terms coupling the left and right nodes of each rail element.
Eq. (1) is solved directly in the time domain using the Newmark method [30] . Application of this method leads to the following time-stepping scheme where the upper index j refers to the time instant t = j t ( t is the constant time step; the presence of the upper index means the variable is time dependant), K is the equivalent dynamic stiffness matrix given by and the constants a 0 to a 5 are given by the following expressions [31] : The constants β and γ are the Newmark parameters, and are set at β = 0 . 25 and γ = 0 . 5 , which leads to unconditionally stable solvers, at least in the case of linear and finite systems. In the current case, the existence of the boundary forces f l and f r , which simulate semi-infinite domains, violates the finite system criteria, meaning that the method does not necessarily have to be unconditionally stable. The choice of these Newmark parameters assumes that within the time interval t the accelerations remain constant and equal to the average between the previous and the next time-step. To solve the time-stepping algorithm in Eq. (2) for u j , ˙ u j and ü j , the interaction forces f l and f r must first be determined. If these forces were resulting from some known excitation, of which one knew their time evolution, then the system could be solved as it is. However, these forces are dependant on the response of the boundary in such a way that the displacements of the three regions must be compatible at the interfaces, and thus these terms must be further expanded. That is addressed in the following subsections. Fig. 4 shows the left and right regions isolated from the mid-region. They are acted upon by the interaction forces f l and f r on the corresponding sides (now acting in the opposite direction: action-reaction pair), and by other external forces f ext , l and f ext , r (not represented in Fig. 4 ), with forces moving along the rail being the most common external excitation. The response of the boundaries of these domains can be written as

Left and right regions
where H α (t) is the impulse response matrix of the left ( α = l ) and right ( α = r ) boundaries (vertical displacement and rotation of the rail and horizontal and vertical translations of lattice masses at the boundaries due to sudden unit-impulse forces applied at t = 0 and at the same boundaries), and u α, ext (t) is the response of these boundaries due to the externally applied forces f ext , l and f ext , r . The minus sign in front of the integral is due to the action-reaction pair at the interface. Also, u α (t) corresponds to u l and u r in Eq. (1) , since there must be compatibility of displacements at the interface.
Due to the discrete nature of the time solver used for the mid-region, the forces f α (t) can only be known at the intervals Thus, in order to evaluate the integral in Eq. (5) , a time evolution has to be assumed for these forces between the discrete time instants. Here, the approach proposed by Bode et al. [32] is followed. That approach assumes that within a time interval the force remains constant and equal to the mean of the force at the beginning and at the end of the time interval, i.e., When compared to other assumptions for the force evolution, such as linear variation, constant and equal to force at beginning, or constant and equal to force at end, the approximation proposed in [32] is the one that showed more robustness in terms of stability of the solver, and thus is the one adopted in this work. Nonetheless, this approximation (or any other) is not consistent with the assumptions made with the Newmark method used to solve the mid-region (in the mid-region, assumptions are made about the acceleration evolution, while at the side regions, assumptions are made about the evolution of the interaction forces), and that can lead to unwanted reflections at the boundaries. These reflections can be minimized by reducing the time step, as discussed later in Section 4 .
Based on the approximation in Eq. (6) and the expressions in Eq. (5) , the boundary responses at the time instant t = j t are given by which can be rewritten as where matrices F j α are calculated with Eq. (8) splits the interface displacements in three components: the displacements induced by the forces acting on the same time interval; the accumulated displacements due to forces acting on the previous time intervals (history term); and the displacements induced by the externally acting forces. That equation provides the framework to study the left and right regions, but expressions for the displacements u j α, ext and for the impulse response matrices H α (t) and F j α still need to be defined. Due to the semi-infinite nature of the side domains, these quantities can be easier calculated semi-analytically in the frequency domain and later on converted to the time domain. Such is explained in Section 3 , but first the coupling of the side regions to the mid-region is explained, so that the whole problem can be solved.

Coupling the side regions to mid-region
Eq. (8) can be rearranged to isolate f j α on the left-hand side, leading to 5 If Eq. (10) is now inserted into the time-stepping algorithm (2), after some rearrangement, namely moving the factor ( F 0 α ) −1 u j α to the left-hand side, one obtains an expanded time-stepping algorithm of the form where vectors f j l and f j r are defined by and where the new equivalent stiffness matrix ˆ K is given by In short, time-stepping algorithms (2) and (11)

Semi-analytical evaluation of interface matrices and vectors
Section 2 explains the main framework of the model proposed in this work. It relies on the knowledge of the impulse response matrices H α (t) and displacements u α, ext (t) , both referring to the interfaces between side regions and mid-region.
Due to the semi-infinite characteristic of the former, the calculation of the mentioned quantities is not straightforward, and it is easier done first semi-analytically in the frequency domain, and then converted to the time domain. The present section explains how this is achieved, starting with the matrices H α (t) .

Calculation of impulse response matrices
The impulse response matrices H α (t) give the response of all degrees of freedom at the interface of the left ( α = l ) or right ( α = r ) region, due to unit-impulse forces applied at any degree of freedom of the same interface (for the calculation of these matrices, the left/right region stands alone; apart from the applied impulse forces, that interface is force/stress-free). Their frequency-domain counterparts, ˜ H α (ω) , are termed dynamic flexibility matrices and contain the harmonic response of all degrees of freedom due to unit harmonic forces applied at any degree of freedom. The rows of the matrices refer to the degree of freedom at which the response is observed, while the columns refer to the degree of freedom where the force is and matrices F j α , as defined in Eq. (9) , can be calculated via The dynamic flexibility matrices ˜ H α can be calculated semi-analytically based on the expressions derived in [29] . In that reference work, the displacements, internal shear and internal moment of the rail, displacements of the sleepers, and displacements of the lattice are given for arbitrary forces acting on the rail or on the lattice, assuming infinite periodicity. To go from these infinite-domain solutions and obtain those of a semi-infinite domain of the same type (to which the dynamic flexibility matrices ˜ H α refer), one must impose that there is no transmission of forces between the left and right semiinfinite domains that compose the infinite domain. This condition is equivalent to creating a semi-infinite domain with a free boundary (as desired), and can be achieved by adding fictitious forces of appropriate amplitude at the part of the infinite domain that is to be removed (i.e., for the calculation of ˜ H l the fictitious forces are applied to the right of the free boundary; for the calculation of ˜ H r , to the left). This procedure is explained next for ˜ H l (for ˜ H r , the procedure is very similar and does not require further explanation).
Consider a periodically infinite domain as represented in Fig. 1 . At each vertical alignment going through the lattice masses, it is possible to apply unit forces F rail and moments M rail at the rail, and horizontal F lat j x and vertical F lat j z unit forces at the lattice masses, and for each of these loads, the expressions from work [29] can be used to calculate the displacements v rail , rotations θ rail , internal shear s rail and moment m rail of the rail at any position, and also the horizontal x lat j and vertical z lat j displacements of the lattice masses at any vertical alignment ( j stands for the depth index of the lattice mass). Apply the mentioned set of forces at alignment i indicated in Fig. 1 (the forces at the rail should be applied at an infinitesimal distance to the left of the alignment), and collect the displacements of the same alignment in matrix U i,i and the displacements of alignment ii (also represented in Fig. 1 ) in matrix U ii,i (first subscript refers to the alignment where the response is calculated; second subscript refers to where the load is applied). For example, matrix U i,i is of the form (the part of the superscript after the comma in each of the matrix entries represents the load location and direction) Matrix U i,i does not correspond to the dynamic flexibility matrix ˜ H l because the expressions from work [29] do not satisfy the free boundary condition needed to simulate a semi-infinite domain, i.e., at alignment i there is transmission of forces from the left to the right side. These forces, here termed connecting forces , correspond to the internal shear and moment of the rail at alignment i , and to the resultant forces at the lattice masses due to the viscoelastic connections between alignments i and ii (see Fig. 1 ). Collect all these connecting forces in a matrix S i,i of the form (the first subscript refers to the alignment where the connecting forces are calculated, while the second subscript refers to the alignment where the forces are applied) where submatrix F i,i , which contains the connecting forces at the lattice, is calculated with In Eq. (18) , matrix ˜ K i,ii is the dynamic stiffness matrix where all the horizontal and diagonal connections between alignments i and ii are assembled, and matrices U lat i,i and U lat ii,i are the submatrices of U i,i and U ii,i , respectively, that collect only the lattice displacements, i.e., after removing the first two rows.
In order to compensate for the non-zero connecting forces S i,i , a second set of forces (termed fictitious forces ) must be applied at alignment ii such that the combination of the connecting forces caused by the loads at alignment i and those caused by the fictitious forces is zero, i.e., ( O is a matrix of zeros) Matrix A fic , which contains the amplitudes of the fictitious forces needed to satisfy the free-boundary condition of the semiinfinite domain, is a square matrix of which each column corresponds to the fictitious force amplitudes for each load scenario at alignment i . It is obtained by rearranging Eq. (19) , and it reads: The dynamic flexibility matrix ˜ H l of the left domain relating forces and displacements at the free surface of the left semi-infinite domain is now obtained by adding the infinite-domain displacements induced by the fictitious forces to those induced by the loads applied at alignment i , i.e., The calculation of matrix ˜ H r follows the same steps as for matrix ˜ H l , with the alignments i and ii swapped, and its final expression is Note that to calculate the left and right flexibility matrices, it may be necessary to consider different material properties if the left and right semi-infinite domains are different. The expressions for their calculation are still Eqs. (21) and (22) , but for each, the quantities must be calculated with the corresponding material properties.

Alternative approach for calculating ˜ H α
Though the method explained in the previous subsection works well and is solely based on the semi-analytical expressions from work [29] , the fact that these expressions require the numerical evaluation of integrals over wavenumbers leads to time consuming calculations. Alternatively, a faster method, which takes advantage of the periodic characteristic of the structure, can be applied. It relies on the discretization of the rail and recursive cloning of cells, leading to some numerical approximation (and propagation of errors). Nevertheless, since the discretization applied is quite fine compared to the characteristic wavelength in the rail (the focus is on low frequencies, below 5 kHz), the introduced errors are minimal and, thus, this alternative approach leads to virtually the same results as the method presented in the previous subsection (which introduces no discretization errors).
Consider the portion of the periodically infinite domain represented in Fig. 1 that is limited by the alignments ii and iii . After discretizing the rail into small beam elements (say, with length l = d/ 2 , where d is the characteristic distance of the lattice -see Fig. 1 ), the frequency-domain equations of motion of the resulting system (lattice masses, sleeper mass and discretized nodes of the beam) are written as where matrices ˜ K αβ ( α, β = l , m , r ) are the dynamic stiffness matrices ( ˜ K αβ = −ω 2 M αβ + i ω C αβ + K αβ ) and where ˜ f ext , α are the external forces applied at each node. For convenience, the degrees of freedom have been divided into those belonging to the left edge (alignment ii , denoted with l ), those belonging to alignment right edge (alignment iii , denoted with r ), and those belonging to the middle (denoted with m ). Taking into consideration that the middle nodes are free of external forces, i.e., ˜ f ext , m = 0 , then the middle degrees of freedom can be removed from the system of equations, obtaining a new relation of the type where the submatrices ˜ K 0 αβ ( α, β = l , r ) are obtained through the expression In Eq. (25) , δ αβ is the Kronecker delta ( δ αβ = 1 if α = β; δ αβ = 0 if α = β), the superscript 0 in matrix ˜ K refers to step zero (at step n , 2 n cells are considered; the recursive procedure is explained next), and the superscript 1 at variables ˜ u and ˜ f refers to the cell number ( 2 0 = 1 ). Matrix ˜ K 0 relates the forces and displacements of the edges of one cell; two of such cells can be placed next to each other and coupled with the portion of the domain limited by alignments i and ii , obtaining in this way a new domain comprising two cells. The motion of this two-cells domain is described by the following system of equations which relates forces and displacements at the extremities of the two cells. Matrix ˜ K couple represents the dynamic stiffness matrix of the alignment coupling the cells, and is basically the same as matrix ˜ K i,ii described in the previous subsection ( Eq. (18) ), augmented with the degrees of freedom (and associated values) corresponding to the nodes of the rail (one on each side). Assuming again that there are no forces applied in the middle alignments, i.e., ˜ f 1 ext , r = 0 and ˜ f 2 ext , l = 0 , then the variables ˜ u 1 r and ˜ u 2 l can be removed from the system of equations, obtaining in this way a new system of the type If this procedure is applied recursively with the newly obtained matrix n times, one obtains a system of equations of the type relating the displacements of the left alignment of the first cell and displacements of the right alignment of the cell 2 n with the forces applied at the same alignments. Thus, by applying this recursive procedure n times, one is able to describe a domain containing 2 n cells with a system of equations containing only the degrees of freedom at the edges of the domain. The expressions to recursively calculate ˜ K n αβ ( α, β = l , r ) are Since each alignment contains few degrees of freedom (for most cases not reaching the hundred figure), each evaluation of (29) is very fast, and within a fraction of seconds thousands of cells can be simulated. Also, as long as there is some damping (in practice, that is always the case), matrices ˜ K n lr and ˜ K n rl will go to zero, and matrices ˜ K n ll and ˜ K n rr will no longer change, meaning that the left and right edges of the modelled domain become uncoupled. In other words, as sensed from each edge, one is effectively modelling a semi-infinite domain. In this way, after convergence of the recursive scheme described in Eq. (29) , the dynamic flexibility matrices ˜ H l and/or ˜ H r are calculated with This recursive procedure revealed to be quite fast, and for the example shown in Section 5 , it takes less than one second to calculate the dynamic flexibility matrices for each frequency, with n ranging between 11 (2048 cells) and 15 (32,768 cells) for convergence to be reached; the other approach takes more than one minute. The relative error between the two approaches was in the order of 10 −4 , meaning that the two procedures will give the same values up to the third significant digit. Now that the dynamic flexibility matrices have been derived, the only unknown components left in Eq. (12) are the displacements u j α, ext at the edges of the semi-infinite regions caused by external loads. These are referred to as the incident displacements and are derived in the following subsection.

Calculation of incident displacements (moving load)
The incident displacements at the edges of the semi-infinite domains are needed to account for the external forces applied at these domains. For the problem tackled in this article, these forces correspond to loads travelling at constant speed, first along the left region and after along the right region. If the incident displacements were not accounted for, then at the moment that the load entered the mid-region, a transient response would be generated, with waves propagating both towards the left and right sides of the interface and polluting the results (the same would happen when the load exited to the right domain). Of course, it can be argued that it is possible to get rid of the transient responses by making the mid-region long enough, but that would defeat the purpose of the approach presented here, which is to reduce the number of degrees of freedom and therefore the simulation time.
Similarly to the procedure employed for the impulse response matrices, the incident displacements u j α, ext are obtained from their frequency counterpart, ˜ u α, ext (ω) , through the expression In turn, the frequency-domain incident displacements ˜ u α, ext (ω) can be calculated by combining the responses induced by a load moving on an infinitely periodic domain with those of point loads, in such a way that the free boundary is respected. While the response of the infinite domain due to a moving load can be evaluated easily, the response due to point loads still requires the evaluation of integrals over wavenumbers, rendering that method less attractive. However, the stiffness matrices ˜ K n ll and ˜ K n rr derived in the previous subsection can be used to obtain the response of the infinite domain to point loads while avoiding the integral evaluation. The whole process is explained next for ˜ u l , ext (ω) . Fig. 1 . Consider that this load crosses the alignment i (for the left domain that is where the free boundary has to be imposed) at time t = t l . This load induces displacements ˜ u ML i at alignment i and displacements ˜ u ML ii at alignment ii (these vectors contain both responses of the rail and of the lattice masses at the corresponding alignments, as explained in Section 3.1 ; expressions for these responses can be found in reference [29] ). Alongside with the displacements, the moving load also induces internal forces at alignment i ( connecting forces between the left and right sides of the alignment), which can be calculated via

Consider a vertical force moving with constant speed V along the rail of the (infinite) domain depicted in
Matrices ˜ K couple αβ are the same as explained in Section 3.2 , and vector ˜ f ML i contains both shear forces and moments at the rail, and connecting forces at the lattice between alignments i and ii , as explained in Section 3.1 . Because matrices ˜ K couple αβ contain stiffness components of a discretized Euler beam representing the rail, some approximation errors will be introduced, but these are very small due to the long wavelength of the rail bending compared to the discretization length. Now, in order to impose the free boundary condition at alignment i , ( fictitious ) forces of amplitudes a fic must be applied at alignment ii such that the combined effect results in the required boundary condition. This step is the same as explained in Section 3.1 , with the difference that now a fic represents a vector of unknown amplitudes and not a matrix, due to the fact that only one loading scenario is considered (the moving load). These fictitious forces induce displacements ˜ u fic i and ˜ u fic ii at alignments i and ii , and connecting forces ˜ f fic i given by As mentioned already, due to its computational cost, it is not desired to make use of expressions in reference [29] to calculate the displacements induced by point forces. Instead, the matrices ˜ K n ll and ˜ K n rr described in the previous subsection are added to the stiffness matrices ˜ K The free boundary condition requires that there is no transmission of forces between left and right sides of the domain at alignment i , which is achieved by the identity ˜ f ML i + ˜ f fic i = 0 . Taking into consideration Eqs. (32) and (35) , that identity is satisfied if For the incident displacements at the right boundary, ˜ u r , ext , the procedure is very similar, with the difference that alignment i is switched with alignment ii , and that the displacements ˜ u ML i and ˜ u ML ii are calculated for the moving load crossing the alignment ii at the time instant t = t r = t l + L m /V , where L m is the length of the mid-region. The final expression for ˜ u r , ext is

Example 1 -Validation
The purpose of the first example shown in this work is to verify all steps of the procedure. For that reason, a simple railway structure, consisting of rails, sleepers and ballast resting on a rigid foundation is considered. Also, in order to compare the results obtained using this method with those of reference [29] , no stiffness variation is considered. The bending stiffness and unit mass of the rail are E I rail = 12 . 08 × 10 6 N . m 2 and m rail = 120 kg (corresponding to two UIC60 rails), the width, mass and rotational inertia of the sleepers are W s = 0 . 3 m , M s = 315 kg and J s = 4 . 73 kg . m 2 , the sleeper spacing (centre to centre) is L = 0 . 6 m , and the thickness of the ballast layer is H ballast = 0 . 3 m . Rail pads are applied between the rail and the sleepers, with vertical and rotational stiffnesses K v = 10 8 N / m and K θ = 2 . 14 × 10 5 N . m (these values corresponds to the pads under both rails, and are based on the values reported in [33] ; the rotational stiffness was obtained assuming that the pad is equally distributed along the contact with the sleeper, of length W p = 0 . 16 m , resulting in K θ = K v W 2 p / 12 ), and corresponding viscous damping C v = 10 5 N . s / m and C θ = 4 . 25 × 10 2 N . m . s . In addition, under-sleeper pads are used between sleepers and ballast, totalling a vertical stiffness of K usp = 1 . 72 × 10 8 N / m per sleeper [34] and corresponding viscous damping C usp = 1 . 72 × 10 5 N . s / m (the track width is assumed to be 2.6 m, same as sleeper length). The damping values of the pads assume a viscous damping ratio of C/K = 0 . 001 s −1 . In the frequency domain, the viscous damping is considered via complex stiffnesses of the type ˜ The properties of the ballast are inferred from the box test results reported by McDowell et al. [5] . There, the authors subjected a ballast sample whose particles sizes ranged between 25 and 50 mm (the median was approximately 30 mm) to successive compressive loads, and observed a ballast stiffness of about 300 kPa/mm. From the reported granulometry, a particle size of d = 0 . 03 m was chosen; this results in 11 lattice masses per column (the bottom one being fixed), N = 11 lattice masses in contact with each sleeper (the stiffness of the under-sleeper pads must be distributed by these masses) and M = 9 free masses in between sleepers. On the other hand, the properties of the particle connections were determined via an inverse problem, in which the springs of a lattice layer with the same thickness as the sample used in the box test (300 mm) were tuned to provide the reported vertical stiffness; together with assuming a Poisson ratio of ν = 0 . 2 [27] and a viscous damping ratio of C/K = 0 . 001 s −1 , the following values were obtained: normal lattice stiffness and damping K n axi = 9 . 68 × 10 7 N / m and C n axi = 9 . 68 × 10 4 N . s / m ; shear lattice stiffness and damping K s axi = 1 . 21 × 10 7 N / m and C s axi = 1 . 21 × 10 4 N . s / m ; diagonal lattice stiffness and damping K n diag = 4 . 23 × 10 7 N / m and C n diag = 4 . 23 × 10 4 N . s / m . Regarding the mass properties of the ballast, an homogenized density of 1800 kg / m 3 has been assumed, which translates into a particle mass of m b = 4 . 21 kg (this may seem a large number, but this value accounts for all the particles at the same level along the width of the track, which is 2.6 m). For the definition of the lattice parameters and rail structure parameters, the reader is referred to the previous work [29] . As can be seen from Figs. 5-6 , the correspondence between blue and red lines in subfigures a-b is to the point that no differences can be detected between the responses obtained with Eq. (15) and with the time domain method. There are nevertheless some small differences between the two, which can be further minimized by further decreasing the time step. The good correspondence between the results obtained with the different methods verifies the validity of this intermediate step.

Verification of impulse response matrices
By comparing Figs. 5 a,c with Figs. 5 b,d, it is observed that the response is larger for the left boundary than for the right boundary. This is because the span between the free edge of the rail and the closest sleeper is larger for the left boundary than for the right boundary (see Fig. 4 ), resulting in a more flexible interface. There are also differences regarding the horizontal response of the lattice mass in Fig. 6 (easier to see when comparing red lines in subfigures c and d), which are justified by the asymmetric cut decided for the boundaries.

Verification of incident displacements
To verify if the procedure to calculate the incident displacements at the boundaries is correct, the results obtained with the process explained in Section 3.3 are compared with those obtained with the same time domain method as in Section 4.1 . The external load is assumed vertical, of unit magnitude and travelling along the rail at constant speed V = 60 m / s . Fig. 7 shows the comparison for the vertical deflection of the rail at the boundaries, and Fig. 8 shows the comparison for the horizontal response of the middle lattice mass.
Once again, the match between blue lines and red lines in Figs. 7 a-b and 8 a-b reveals a very good correspondence between the method described in Section 3.3 and the time domain method, thus verifying the correctness of the former. In addition, by comparing the responses of the left boundary and those of the right boundary, very distinct behaviours can be discerned that are not justified solely by the asymmetry of the left and right domains. The reason for these big differences is in the moving nature of the load; while in the left domain the load is moving towards the free edge and then leaves it, in the right domain the load suddenly enters it and then moves away from the free edge. That is why the response of the left boundary increases gradually and then shows a free decay behaviour (starting the moment the load crosses the free edge), and why the response of the right boundary is zero until the load reaches the edge, then presents a sudden increase, and decays afterwards.

Verification of the whole procedure
The whole process as explained in Section 2 is verified next. With that intention, the responses of an infinite periodic track due to a load moving at the constant speed of V = 60 m / s and unit magnitude obtained using the proposed method and using the equations presented in work [29] are compared. For the proposed method, the mid-region consists of ten identical cells, and for the time step, two cases are considered: t = 0 . 001 s (as in the previous examples) and t = 0 . 0 0 01 s . The results to be compared are the vertical deflections of the rail at the contact points with the sleepers. Because the magnitude of the moving load is constant in time, the response of any given point at two distinct cells must be the same (apart from a time shift), and thus, the expressions in [29] are used to calculate the response of the rail only at one position. Fig. 9 shows the results obtained with the two methods; Figs. 9 a and 9 c show the responses at the different positions over time, as obtained with the proposed method, for the time steps t = 0 . 001 s and t = 0 . 0 0 01 s , respectively, and Figs. 9 b and 9 d show the same responses for the same time steps, but with time shifts such that the periodicity of the response can be assessed. Figs. 9 b and 9 d also depict in dashed black line the theoretical response as obtained with the expressions from [29] . The vertical dashed lines in Figs. 9 a and 9 c represent the instants at which the moving load enters and leaves the mid-region.
From Fig. 9 a it is observed that at the moments that the load enters and leaves the mid-region, some disturbances occur that affect the results. This is also observable in Fig. 9 b, where it can be seen that the response of the rail at the sleeper positions is not perfectly periodic, i.e., there are some disturbances around the expected response (in dashed black line). The reason for these disturbances is the incompatibility between the assumptions made for the left and right regions and those made for the mid-region; while for the former it is assumed constant force during one time interval, for the latter it is assumed constant acceleration. By reducing the time step, it is expected that the inconsistencies become less relevant, and such is proven in Figs. 9 c-d, where it is seen that no disturbances can be seen when the load enters or leaves the mid-region, and also that there is a virtually perfect overlap of the responses of the rail at different positions. Despite the differences between expected and obtained responses for t = 0 . 001 s , the maximum deflection is well captured (only at the first sleeper the maximum deflection is slightly smaller) and the overall shape as well. Thus, this time step is used in the example solved in the next section, so that the calculations do not become too time consuming.

Example 2 -culvert passage
To demonstrate the capabilities of the model, the case of a railway line passing over a culvert is considered next. It represents a transition from a soft foundation to a very stiff foundation and again to a soft foundation, and is based on an existing stretch of the Dutch railway line near the city of Gouda, which is thoroughly described in previous publications by other authors [35][36][37] . Here, only the parameters that are essential for the simulation are described, and the modelling assumptions are explained. Fig. 10 schematizes the conditions at the site: a culvert, with longitudinal dimension corresponding to four sleepers, is placed at the depth of 60 cm below the ballast layer; the culvert rests on piles, making it practically immovable (rigid); two concrete slabs with thickness of 30 cm and length corresponding to 7 sleepers lie on each side of the culvert with the intention of mitigating the abrupt stiffness variation; the foundation consists of a sand embankment on top of clay and sand layers.
For the superstructure (rails/sleepers/ballast), the parameters and modelling strategies (lattice with characteristic diameter d = 3 cm ) are the same as explained in Section 4 . Regarding the underlying foundation, it is modelled as a lattice with the same characteristic diameter d as the ballast and whose properties are discussed in the next subsection. For the culvert, since it is practically immovable, it is modelled by constraining the motion of all lattice masses that are at the location of the culvert (thus simulating a rigid and immovable body). In turn, the approach slabs are simulated by assigning the properties of concrete to the lattice masses and lattice connections at corresponding locations. The properties assumed for concrete are Young's modulus E = 20 GPa , Poisson's ratio ν = 0 . 2 and unit density ρ = 2300 kg / m 3 , and based on these values and assuming that the slabs are 2.6 m long in the transverse direction (same as track width), the following lattice properties are   [27] : particle mass m = 5 . 38 kg , normal spring K n axi = 3 . 98 × 10 10 N / m , shear spring K s axi = 3 . 61 × 10 9 N / m , and diagonal spring K n diag = 1 . 81 × 10 10 N / m . It is assumed that the slabs are hinged to the top of the culvert.

Foundation properties: from 3D to 2D
The foundation consists of a sand embankment on top of a layered soil. The mechanical properties of embankment and underlying soil are indicated in Table 1 . The properties in Table 1 were determined assuming real conditions, i.e., assuming 3D propagation of waves. However, the model developed in this work is 2D, meaning that the propagation of waves along the transverse direction is disregarded. The 2D assumption is reasonable for the superstructure, since there most of the energy propagates longitudinally or vertically, but less valid for the foundation, where energy propagates also along the transverse direction. In this way, one cannot simply assume plane-strain conditions for the foundation and use the parameters in Table 1 to determine the corresponding lattice properties, otherwise the results would be affected by the lower decay rate of waves/energy. Instead, one must find equivalent parameters for a 2D model which lead to a similar response as that obtained under 3D conditions. The procedure used to derive the equivalent properties is explained next.  The ballast layer rests on the foundation along a certain width on which the interaction forces are non-uniformly distributed. These interaction forces lead to displacements at the interaction surface, which are also not constant. In order to condense all this information into a plane (remove the transverse direction), certain assumptions must be made about the distribution of the interaction forces and the displacements at the interface. In this work, it is assumed that the forces are uniformly distributed along the width of the interface, and that the motion at the interface can be described by the displacement at the middle of the interface (other assumptions could have been made [38] ). Based on these assumptions, it is possible tune the foundation properties of a plane-strain model such that the force-displacement relationship at the surface matches (as much as possible) the same relation of the 3D model.
For the present example, the tuning of the 2D model is based on the force-displacement relationship induced by loads moving at 50, 100, 150 and 200 km/h (i.e., the steady-state response of the plane-strain model will be matched, as good as possible, with that of the 3D model). In addition, in order to save computational time, the 2D model consists of only one layer that is 1.2 m thick and rests on a viscoelastic foundation ( Fig. 11 depicts the 3D scenario and the equivalent 2D scenario to be tuned). The single layer is assigned the shear wave speed C S = 140 m / s (just like the sand embankment) and the Poisson's ratio ν = 0 . 25 (not the same as the sand embankment because the lattice is limited to Poisson's ratio of 0.25 [39] ; above that value, the stiffness of the shear spring of the lattice becomes negative, which leads to a divergent response), and the remaining parameters are obtained from the tuning process, for which the Matlab function fminsearch is used to minimize the integral of the absolute differences between the 3D and 2D responses. The following tuned parameters are obtained: density of the layer ρ = 3560 kg / m 3 ; viscous damping ratio η = C/K = 8 . 35 × 10 −4 s −1 (because the lattice model is solved in the time domain, damping is considered viscous instead of hysteretic, as defined in Table 1  As can be seen in Fig. 12 , the match between the results of the tuned 2D model and those of the 3D model is not perfect; the displacement under the load is not far from the desired value, but the attenuation rate of the response is higher for the 2D case than for the 3D case. A better match can be obtained if more parameters are used in the tuning process (for example, by adding more layers whose properties can be tuned), but that leads to more time consuming calculations due to the added degrees of freedom (not to mention the added difficulty in the tuning process). Also, matching the response for moving loads as seen at the surface does not necessarily mean that the response due to point loads at varying frequency or at other depths will also show good correspondence. Thus, it may be a more robust procedure to, instead, match the transfer functions in the wavenumber-frequency domain at multiple depths, not just the surface. Nevertheless, for the purpose of this work, which is to introduce the method discussed in Sections 2 and 3 and to demonstrate its applicability, the values obtained in the previous paragraph are accepted. From those values, the following parameters for the lattice representing the viscoelastic layer are calculated [27] : particle mass m = 3 . 20 kg , normal spring K n axi = 1 . 40 × 10 8 N / m and damping C n axi = 1 . 17 × 10 5 Ns / m , shear spring K s axi = 0 and damping C s axi = 0 , and diagonal spring K n diag = 6 . 98 × 10 7 N / m and damping C n diag = 5 . 83 × 10 4 Ns / m . The values of the viscoelastic foundation must be lumped at each lattice mass at the lower boundary, which is achieved by multiplying the vertical and horizontal stiffness and damping values by d.

Results
The external force used in this example corresponds to a point load of unit magnitude moving along the rail at the constant speed of V = 60 m / s (220 km/h). Three scenarios are considered for comparison: the case in which the culvert is not present; the case in which the culvert is present but no approach slabs are considered; the scenario depicted in Fig. 10 , in which both culvert and approach slabs are considered. The total number of cells (sleepers) modelled in the mid-region amounts to a total of 28, of which 10 correspond to the regions right before and after the approach slabs (5 on each side; in principle these 5-cell buffers are not necessary, they are added for result visualization), 14 correspond to the approach slabs (7 for each slab), and 4 correspond to the culvert. The responses to be compared are the deflections of the rail at the contact point with the sleepers ( Fig. 13 ), the maximum force transmitted by each sleeper to the ballast ( Fig. 14 a), the energy dissipated by the ballast layer under each sleeper ( Fig. 14 b) and the energy transmitted to the ballast layer by each sleeper ( Fig. 14 c). The energy dissipated by the ballast under a given sleeper is calculated by adding the energy dissipated by all the ballast viscoelastic connections between two consecutive mid-span positions. Note that the blue lines in Fig. 14 refer to the scenario with no stiffness variation and therefore should be straight and horizontal (each cell behaves the same); the oscillations of these lines are caused by the differences in the approximations assumed for each domain (as discussed in Sections 2 and 4 ), and can be reduced by decreasing the time step. Fig. 13 a shows that in the absence of approach slabs and culvert the response of the rail is the same for every cell (apart from the small differences already discussed in Section 4 ), a fact that is expected due to the periodic nature of the structure. On the other hand, Fig. 13 b shows that the presence of the culvert leads to a reduction of the rail deflection around its location, which is also expected due to the increase of the support stiffness caused by the culvert. The reduction of the maximum rail deflection is less abrupt in the presence of the approach slabs, as seen in Fig. 13 c, suggesting that the approach slabs indeed provide a gradual transition from the soft to stiff foundations. Interestingly, in the latter two cases (with and without approach slabs), it only takes a couple of sleepers before and after the culvert or the approach slabs for the deflection of the rail to go back to the value of the steady state.
The maximum force transmitted by a sleeper can provide an indication of how much settlement might occur under that sleeper; the higher the transmitted force, the larger the settlement [41] . In that respect, Fig. 14 a suggests that the maximum transmitted force increases in the presence of the culvert (red line); at the sleepers right before and after the culvert the force decreases by about 1.5%, but above the culvert the force increases by about 6%. The presence of the approach slab reduces the increase to about 5% (black line), but the increase is spread over a larger distance (starting and finishing at the edges of the approach slabs), meaning that settlements may occur over a wider area of the track.
On the other hand, some authors suggest that the damage in track components is dependant on the energy dissipated by the component [ 12 , 42-44 ], and based on that the energy dissipation plots in Fig. 14 b can give an indication of the sections of the track more susceptible to ballast settlement. The figure indicates that less energy is dissipated at the sections with increased stiffness (culvert and approach slabs) than at the homogeneous sections of the track (the decrease in energy dissipation is in the order of 20%). This suggests that less ballast settlement is to be expected at these regions of larger stiffness, which is conflicting with the observation of the previous paragraph.
The energy flow from the sleepers to the ballast (from the upper structure to the foundation) does not give a direct indication of the potential settlement, since once transmitted to ballast, energy does not necessarily have to dissipate there, as it can be further transmitted to the foundation or propagate horizontally (radiation). However, one can try finding a correlation between energy transmitted and energy dissipated to make indirect assessments of settlement. If such correlation exists, then one can simplify the models by condensing the ballast/foundation layers into lumped connections, and from these simplified models assess the damage based on energy transmitted. Unfortunately, the comparison of Figs. 14 b and 14 c does not reveal any relation between the two, and so not much more can be inferred from Fig. 14 c other than that the presence of culvert and approach slabs significantly alter the energy flow within the system; while for a homogeneous track the energy always propagates downwards (positive values), for a track with stiffness variation the energy can also propagate upwards (negative values, before the transition) as a result of reflections at the stiffer domains.
The results observed in the current section are valid only for the properties and load speed considered. Also, they do not account for the interaction between the train and the track, since the moving force is constant in magnitude. The inclusion of a mass-spring oscillator representing an axle/bogie system will lead to time varying force amplitudes, which in principle are amplified around the region of stiffness variations, and thus might result in force distributions different from those reported in Fig. 14 a. Also, no damage mechanism has been included in these simulations, and thus the energy dissipation plots presented in Fig. 14 b must be taken lightly, since non-linear mechanisms such as frictional sliding will alter how the energy will propagate and how much energy will be dissipated in a given connection. Despite the current limitations of the model, and the argument that the results reported in this section to some extent lack a realistic loading scenario, the contradictory inferences from Figs. 14 a and 14 b regarding ballast settlement lead to the conclusion that at least one of the linear indicators used in this work (maximum force and dissipated energy; termed linear because they were obtained from linear analyses) does not provide correct assessments regarding ballast settlement. Once a non-linear description of the lattice and the interaction with a moving oscillator are included, more reliable results can be obtained from which proper settlement/damage indicators may be derived.

Concluding remarks
This paper presents a method for analysing the response of ballasted railway tracks containing stiffness variations in the support, with emphasis on modelling the ballast with a lattice layer, whose dispersion characteristics differ from those of a continuum. The main advantage of the proposed method is that, unlike other FEM based methods, the boundaries are nonreflective and therefore can be placed close to the stiffness variation without polluting the results with unwanted reflections. This means that the model can be limited to the region of interest, which in the context of ballast settlement is the zone of stiffness variation plus a few sleepers before and after the transition.
As presented, a given limitation of the model is that the Poisson's ratio of the foundation (lattice) is limited to 0.25. Above that value, the stiffness of the shear springs of the lattice become negative, which is non-physical and leads to divergence of the solution. Though that might not be an issue for the ballast, since its Poisson's ratio is about 0.2, for the foundation, and namely for saturated clays or sands, that will be an issue, since their Poisson's ratios are near 0.5 and 0.3, respectively. A way around this issue can be coupling the lattice layer representing the ballast with a continuous medium representing the underlying foundation, as done by Ricci and collaborators [45] . That approach will make it more difficult to compute semi-analytically the incident displacements and the matrices of the transmitting boundaries as described in Section 3 , but if combined with periodic FEM models, can become a very powerful way around the issue.
More attention is needed when defining the properties of the foundation. Since the model is 2D, intrinsically, it cannot capture the spread of energy along the transverse direction, and thus equivalent 2D properties must be defined. The method used in this paper to define the equivalent properties was a simple one, in which the 2D displacement field induced by moving loads was matched to the corresponding 3D field, but better strategies may already have been reported in the literature. Because the purpose of this article is to introduce the model, focus was not put on the exact definition of the parameters, but that does not interfere with the applicability of the method.
As a final remark, since the method is formulated in the time domain, the model can be easily adjusted to incorporate nonlinearities (ballast settlement, hanging sleepers, etc.). This will be addressed in a future publication, in which a constitutive model for the lattice representing the non-recoverable behaviour of ballast will be presented. In its current state, the model can be useful to understand instantaneous dynamic behaviour of the track in a transition zone, and analyse features like force amplification, energy propagation/dissipation and vehicle-track interaction. After augmented with these non-linear capabilities and interaction with an oscillator, the model can become a powerful tool for finding proper damage/settlement indicators and for studying the efficiency of mitigation measures.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.