Fast simulation of railway bridge dynamics accounting for soil–structure interaction

A novel numerical methodology is presented to solve the dynamic response of railway bridges under the passage of running trains, considering soil–structure interaction. It is advantageous compared to alternative approaches because it permits, (i) consideration of complex geometries for the bridge and foundations, (ii) simulation of stratified soils, and, (iii) solving the train-bridge dynamic problem at minimal computational cost. The approach uses sub-structuring to split the problem into two coupled interaction problems: the soil–foundation, and the soil–foundation–bridge systems. In the former, the foundation and surrounding soil are discretized with Finite Elements (FE), and padded with Perfectly Match Layers to avoid boundary reflections. Considering this domain, the equivalent frequency dependent dynamic stiffness and damping characteristics of the soil–foundation system are computed. For the second sub-system, the dynamic response of the structure under railway traffic is computed using a FE model with spring and dashpot elements at the support locations, which have the equivalent properties determined using the first sub-system. This soil–foundation–bridge model is solved using complex modal superposition, considering the equivalent dynamic stiffness and damping of the soil–foundation corresponding to each natural frequency. The proposed approach is then validated using both experimental measurements and an alternative Finite Element–Boundary Element (FE–BE) methodology. A strong match is found and the results discussed.


Introduction
The response of railway bridges is strongly affected by soil-structure interaction (SSI), especially under resonance conditions (Romero et al. 2013). Radiation and material soil damping influence the modal parameters and usually mitigate the structural vibration levels. However, the interaction between the sub-structure and the soil is seldom considered in numerical models when solving the dynamic problem. In order to predict realistically the vibrational response of bridges under railway traffic, or to asses the structural integrity for new operational situations, SSI should be considered. In the design of new structures this may lead to optimized alternatives and in the case of existing ones, it will permit to assess the bridge performance when facing operational challenges (e.g. increase in the speed or capacity of services) and improve model calibration procedures.
Several researchers have investigated the dynamic response of the bridge accounting for SSI. The soil surrounding the bridge foundations can be modelled applying the Boundary Element method (Domínguez 1993) or the FE method adding non-reflecting boundaries (Lysmer and Drake 1972). Galvín and Domínguez (2007) studied the vibrations induced by high-speed trains crossing a curvlet-type railroad underpass using a three-dimensional (3D) BE model. Báez et al. (2018) analysed the aforementioned structures following a substructuring approach, and Zangeneh et al. (2018) by means of a FE model. In all cases, the authors found the influence of SSI relevant. Romero et al. (2012Romero et al. ( , 2013 studied the effect of SSI on the resonant response of railway bridges. They concluded that the bridge fundamental periods and damping ratios increase when soil interaction is taken into account. Later, Doménech et al. (2016) and Martínez-Rodrigo et al. (2018) analysed SSI effects over an extensive catalogue of bridge typologies prone to experience high transverse accelerations at resonance for modern trains and high design velocities. They studied the structural response using a coupled 3D BE-FE model and concluded that incorporating SSI effects in the analysis of new simply supported (SS) bridges or in the analysis of existing ones for retrofitting purposes could be of great importance in the assessment of the Serviceability Limit Sate (SLS) associated to traffic safety. González et al. (2020) studied the performance of inclined pile foundations on the seismic response of bridges. The authors computed impedance functions and kinematic interaction factors of the pile foundations and compared them with the results of two different approaches: a Winkler-type model and a BE-FE model. Zangeneh et al. (2019) presented a simplified discrete model for calculating the modal parameters of the fundamental bending mode of railway beam bridges on viscoelastic supports. They proposed exact, closed-form expressions for the fundamental frequency and modal damping ratio, considering the effect of the dynamic stiffness and dissipation capacity of the foundation-soil system. Recently,  investigated the effect of the foundation mass on the fundamental modal properties of the beam bridges. Gara et al. (2019) studied the effect of the foundation-soil system on the interpretation of experimental ambient vibration tests and on the calibration of numerical models. They concluded that the common practice of updating fixed base numerical models to fit experimental results should be prudently evaluated in the case of bridges where SSI is significant.
In the aforementioned research, simplified foundation geometries are mostly considered due to the high computational cost of the chosen numerical approaches. However, the computational effort could be reduced to a certain extent during SSI simulations using Perfectly Matched Layers (PML) to limit the size of the soil domain (Basu andChopra 2003a, b, 2004;Basu 2009;Kausel and de Oliveira Barbosa 2012). For example, Davoodi et al. (2018) studied the response of simplified geometries of different foundations using the absorbing layers provided by the PML incorporated in a FE model. Lopes et al. (2014) studied the problem of a track-tunnel-ground system using a two-and-a-half-dimensional (2.5D) FE-PML approach. Esmaeilzadeh Seylabi et al. (2017) solved the wave propagation problem in the frequency-domain on a FE model applying PML to numerically simulate measured dynamic impedances of soil-pile systems.  investigated the effect of dynamic SSI on the resonant response of portal frame bridges founded on a rigid base, and the parameters governing their interaction. The authors compare the results given by a simplified model with those from a full 3D FE-PML model. Figure 1 shows a schematic representation of the railway beam-bridge coupled problem. The problem is usually divided into two sub-domains. First the bridge deck, supports and foundations are commonly represented using FEs. Then this structural domain is coupled to the soil subspace. Alternatively, if the complete problem is modelled, and therefore modal superposition is not applied, the corresponding dynamic equilibrium equations can be solved using either implicit or explicit integration schemes, in the time or frequency domains. These approaches require large amounts of computational resources and time since, in addition, it is usually required to obtain the bridge response under several train types and speeds. An alternative and faster approach is to decouple the problem and to compute the structural response as the superposition of the first modal contributions. However, this transformation is not straightforward when the soil is accounted for, due to the presence of non-proportional damping.
In this paper a general approach is proposed to solve the SSI problem in railway bridges applying modal superposition. A sub-structuring approach involving PML permits the modelling of arbitrary foundation geometries. First, the problem in Fig. 1 is decoupled into (i) the soil-foundation problem (Fig. 2), and, (ii) the soil-foundation-bridge interaction problem (Fig. 3). In the first, the equivalent frequency dependent dynamic stiffness and damping characteristics of the soil-foundation system are computed using the FE-PML method. In the second, the dynamic response of the structure under railway traffic is computed using a FE model that includes spring/dashpot elements at the supports, where the equivalent properties at the bridge-foundation connection points are extracted from the first problem. This soil-foundation-bridge model is solved by modal superposition considering the equivalent dynamic stiffness and damping of the soil-foundation corresponding to each natural frequency. To do so, a model updating procedure is implemented. As the model presents non-proportional damping, the complex eigenvalue problem is solved in order to compute the natural frequencies and mode shapes.
The study considers the vertical force transmission between the deck and the foundation, as the aim of the investigation is to develop and present the methodology. Nevertheless, the approach could easily be extended to more complex bridge-foundation-soil systems.
The outline of the paper is as follows. In Sect. 2, the numerical formulation is presented. In Sect. 3 the modal properties and dynamic response of a particular soil-foundation-bridge system are validated and compared to those from previous publications which Soil Foundation Bridge P Fig. 1 Schematic representation of soil-foundation-bridge coupled system are computed applying different numerical techniques. In Sect. 4 the response of a real bridge is reproduced and compared to experimental measurements. Finally, in Sect. 5 the main conclusions are summarized.

Numerical formulation
The soil-foundation-bridge interaction problem ( Fig. 1) is decoupled (Figs. 2, 3) with the aim of implementing a fast and accurate analysis tool that can be used in practice for the design and performance evaluation of these types of structures.

Soil-foundation interaction problem
The dynamic stiffness of the soil-foundation system is computed applying the FE method. In order to represent the absorbing boundaries for the wave equation, PML is used. Following, the formulation of the PML is briefly summarized (Basu and Chopra 2003b;Bindel and Govindjee 2005;ANSYS, Inc. 2021).
The PML layers absorb and attenuate the arriving waves. In the frequency domain, the PML medium is formulated assuming a harmonic time dependence of the displacement, stress where ′ ij are the components of the stress tensor in the local coordinate system, and is the density of the medium.
The components of the strain tensor in the local coordinate system are: The stretching function is where the functions f e j and f p j attenuate evanescent and propagating waves, respectively. If evanescent waves are not present, the stretching function results The attenuation functions used in this work and referred to in Fig. 4 are: The parameters f 0j and L Pj are set to obtain a normal reflection coefficient of −60 dB , i.e.:

Fig. 4 Attenuation function
The versatility of the FE methodology allows the simulation of different foundation types and soil stratigraphies. Firstly the soil-foundation system shown in Fig. 2 is considered in order to compare the proposed methodology with the results presented in Doménech et al. (2016) and Martínez-Rodrigo et al. (2018).
Three homogeneous soil types are defined with flexibilities covering the AASHTO classification (AASHTO 2012). In particular s and p-wave speeds of c s = {150, 220, 365} m/s and c p = 2c s are admitted, meaning that Poisson's ratio s = 1∕3 . Soil density and material damping are s = 1800 kg/m 3 and s = 0.05 , respectively, for all cases. Identical 2B × 2B rigid foundation plates with B = 2.5 m are considered in all the cases initially.
The soil-foundation system is analysed using the FE model shown in Fig. 2 implemented in ANSYS. In this figure only volumes are indicated for simplicity, and not the actual mesh. A massless foundation plate is considered for a better comparison with the cited works, which is modelled with four-node shell elements with six degrees of freedom (dof) per node. The plate is surrounded by a volume of soil represented with eight-node solid elements with three dof per node. At the boundaries of the soil elements, additional solid elements are created incorporating the PML formulation in order to act as absorbing boundaries.
The number of PML layers determines the absorbing efficiency of the PML region. The absorbing efficiency depends on the excitation frequency, the distance between the source and the PML domain, the PML domain thickness L Pj , and the wave propagation velocity. The number of layers of the PML region is obtained from the PML thickness and the required element size. It is recommended a number of layers higher than 2 to reproduce the attenuation distribution. Moreover, it is advisable to model the near-field soil, at least, by two elements (ANSYS, Inc. 2021).
A total soil volume of 20 m × 20 m × 10 m is generated, where the thickness of the PML boundary equals 5 m , according to the mentioned requeriments. Therefore, the distance between the PML region and the foundation is 2.5 m . The size of the elements is selected based on the frequency. In particular, the mesh is discretized in 10 elements per wave length = 2 c s ∕ (Coronado and Gidwani 2016) with a minimum element size of 0.5 m.
When a harmonic vertical force F s ( ) is applied at the center of the foundation plate with frequency , the dynamic stiffness is defined as Domínguez (1993): where v s is the vertical displacement of the foundation. K d v ( ) can be expressed as: The real and imaginary parts of the previous equations are related to the stiffness and inertia, and to the damping properties of the soil, respectively. Therefore, where K v ( ) = Re{K d v ( )} stands for the stiffness, and C v = Im{K d v ( )}∕ represents the damping of the system in terms of the frequency.
(4) 20 log 10 = 20 log 10 10 −3 = −60 Figure 5 shows the vertical compliance for the surface square foundation. The results are normalised to a 0 = B∕c s and to the compliance F 0 = Bc 2 s s ∕K d v as in Domínguez (1993). The normalised results for the three different soils types with c s = {150, 220, 365} m/s match as indicated in Domínguez (1993). The agreement with the results computed with a BE model, shown in the same figure in solid black trace, is good. Figure 6 presents the soil-foundation system stiffness and damping for the same soils. Moreover, the results presented in Gazetas (1991) are superimposed with circles. Gazetas (1991) provided results in the dimensionless frequency range a 0 from 0 to 2, i. e., 0-19 Hz , 0-28 Hz and 0-46 Hz , respectively, for the three soils considered. In the cited research, the stiffness decreases approximately 40% percent and the radiation damping remains fairly constant while the material damping decreases with the frequency in the considered range. The agreement between the proposed methodology and the results presented in Gazetas (1991) is good. At low frequencies, the FE-PML model predicts slightly higher values for the stiffness, and lower for the damping for the two stiffest soils. As the frequency increases, both properties decrease. The decrement of the material damping is higher in Gazetas (1991).
The advantage of the methodology presented herein is that it permits the representation and analysis of any type of foundation (abutments, piles, etc.) and soil stratigraphy. This study focuses on the numerical validation of the approach based on previous results. Cases dealing with more complex geometries will be presented in future works. Figure 3 shows the soil-foundation-bridge decoupled model. Beam bridges of length L, cross-section bending stiffness EI z , mass per unit length m b , and structural damping are discretized using two-node Bernoulli-Euler beam elements, well suited for the analysis of railway bridges in this study given the slenderness ratios and the frequency range of interest from 0 to 30 Hz (Doménech et al. 2016). The spring/dashpot elements at the bridge-foundation connection points, with equivalent stiffness and damping values K v and C v , are obtained from the FE-PML model as explained in Sect. 2.1. Therefore, the crossinteraction between the two foundations of the bridge is neglected. The bridge FE model with N dof is solved using complex modal superposition (Hurty and Rubinstein 1965). The equilibrium equation of the system is:

Soil-foundation-bridge interaction problem
where M , C , and K are the mass, damping and stiffness matrices, respectively. Due to the SSI effect, represented by the spring/dashpot elements at the supports, damping in the problem is non-proportional, i.e., (M −1 C)(M −1 K) ≠ (M −1 K)(M −1 C) . In the case of nonproportional damping, the position of each dof in each mode shape is defined by the amplitude and the phase, which requires 2N equations to determine the solution of the N dof structure. Equation 8 can be rewritten as a system of 2N equations to be solved in this case as: Equation 9 is expressed as a first order matrix equation: where In the previous equations, A and B are real symmetric matrices with dimensions 2N by 2N, and y is the state vector.

Free vibration response
In the free vibration case, Eq. 10 takes the form The non trivial solution of this linear system of equations may be expressed as

Ȧy(t) + By(t) = 0
where s j represents the j-th element of a set of 2N eigenvalues, and j are the corresponding eigenvectors: The natural frequencies, damped natural frequencies and the modal dampings may be obtained from the eigenvalues s j as j = |s j | , dj = |Im{s j }| and j = −Re{s j }∕|s j | . Therefore, in the problem analysed, Re{s j } should be negative and Im{s j } ≠ 0 . Moreover, orthogonality conditions T j A k = 0 and T j B k = 0 for j ≠ k apply as proven in Hurty and Rubinstein (1965), where the superscript T indicates transpose.
The natural frequencies of the soil-foundation-bridge system are computed as per Eq. 12. First, the values of K v and C v obtained from Fig. 6 for the j-th natural frequency of the SS beam are selected as an initial solution. Then, an iterative problem is solved by updating the values of K v and C v for the damped natural frequencies of the bridge, upto the relative error in the n-th iteration, (f n j − f n−1 j )∕f n−1 j < 0.01 . The convergence of the problem is very efficient, involving generally only a few iterations.

Forced vibration response of non-proportional damped structures
The solution to Eq. 10 can be expressed as: Applying orthogonality conditions and normalising the eigenvectors to the matrix A , i. e., T j A j = 1 , Eq. 10 can be uncoupled in a set of 2N equations as (Hurty and Rubinstein 1965) where j = T j B j = −s j and p j (t) = T j P(t). Taking into account that the eigenvalues and eigenvectors for the structures analysed are pairs of complex conjugates, the displacements in the bridge can be computed from Equation 16 may be solved numerically. It corresponds to a non-stiff differential equation that, in what follows, is solved using a Runge-Kuta (4,5) explicit algorithm (Dormand and Prince 1980;Shampine and Reichelt 1997).
The effects of SSI on the seismic structural response are large and they should be accounted for. Carbonari et al. (2011) used lumped parameter models to approximate the frequency-dependent behaviour of the SSI in multi-span bridges and to evaluate their seismic responses. Poul and Zerva (2019) implemented the PML for modeling the unbounded domain in the seismic response evaluation of a concrete gravity dam. Chen et al. (2021) computed the frequency response of offshore monopile foundations to seismic excitation (13) y(t) = j e s j t using the PML as an absorbing boundary condition. The present methodology can be extended to seismic problems applying the corresponding accelerograms ( ü g (t) ) in the force term of Eq. 9, by an additional term − ü g (t) . Then, u(t) = u t (t) − u g (t) , where u t (t) represents the total displacement of the structure. The proposed approach allows the computing of the displacements and the stresses in the bridge.

Numerical analysis
Several of the cases presented by Doménech et al. (2016) and Martínez-Rodrigo et al. (2018) are undertaken following the approach presented in Sect. 2, with the aim of performing a numerical validation and comparing the results with those obtained by applying a different numerical procedure. In the aforementioned works, the authors analyse SSI in beam bridges with the same superficial rigid plate foundations using a fully coupled 3D BE-FE model integrated in the time domain. In this case, the Green's function for an elastic half-space is used as the fundamental solution for the soil, and material damping is approximated as a constant value. In order to compare with the results presented in Doménech et al. (2016), the same structure studied in the mentioned reference consists of a beam bridge of length L = 17.5 m , mass per unit length m = 17,500 kg/m , and fundamental frequency equal to 35% of the simplified method frequency band in Eurocode 1, is analysed under the circulation of the HSLM trains from these standards. A structural Rayleigh damping of 1.18% is assigned according to the design recommendations of Eurocode 1 for prestressed concrete bridges. The three soil types presented in Sect. 2.1 are studied. 5.00 % soil material damping is used in all cases.
A modal analysis of the bridge is performed as per Eq. 12. Non-oscillatory modes with real eigenvalues are excluded from the analysis. The first vertical bending frequency of the bridge in the simply supported (SS) case is f SS 1 = 6.87 Hz . When the interaction with the soil is taken into account, the fundamental frequencies resulting from the proposed approach for the three soil types are: f 150 1 = 6.13 Hz , f 220 1 = 6.43 Hz , and f 365 1 = 6.69 Hz , where the superindices indicate the soil shear wave velocity. Moreover, the modal damping ratios, which include the effect of both material and radiation damping, in addition to structural damping, computed from the complex eigenvalue problem are: 150 1 = 9.44% , = 7.14% . This mode is highly damped, mainly for the softer soil. The soil-foundation system stiffness and damping obtained from the iterative process are shown in Fig. 7 for these particular cases.
It should be also remarked that the second bending natural frequency for the system, accounting for these damping and stiffness values, decreases as the shear wave velocity increases. It can be explained by the parametric analysis shown in Fig. 8. The frequencies corresponding to modes 1 and 2 for the stiffness values shown in Fig. 7a and different C v ∕K v ratios are plotted in solid trace. The frequency values obtained for the C v values obtained from the iterative procedure and presented in Fig. 7b are superimposed with circles. Figure 9 shows the mode shapes of the soil-bridge systems corresponding to the first and second natural frequencies. The real and imaginary parts are plotted, and in the complex mode shapes, each degree of freedom has a different phase angle. Figure 10 shows the free vibration response of the beam at mid-span after the passage of a single load in terms of the non-dimensional speed K SS 1 = v∕(2f SS 1 L) for different soil properties. These free vibration envelopes permit the identification of velocities leading to either local maxima or cancellation conditions (minima of free vibrations). The amplification of the response to be expected at resonance, depends on the amplitude of these free vibrations, which will add in phase in the resonant case. Figure 11 presents the maximum bridge acceleration at mid-span under the loading of HSLM-A1 and HSLM-A6 trains (Eurocode 1 2002). The characteristic distances of these trains are d = 18 m and d = 23 m , respectively. HSLM-A1 and HSLM-A6 theoretically induce a second resonance of the fundamental mode of the bridge at v 1,2 = f SS 1 × d∕2 = 6.87 × 18∕2 × 3.6 = 222.59 km/h , and v 1,2 = f SS 1 × d∕2 = 6.87 × 23∕2 × 3.6 = 284.42 km/h , respectively (non-dimensional speed 0.26 and 0.33). The dashed horizontal line indicates the SLS of the vertical acceleration for ballasted track bridges ( 0.35g = 3.43 m/s 2 ).  In Fig. 11a, b the response of the bridge under the two aforementioned trains computed following the present approach is compared with the results from Doménech et al. (2016). The HSLM-A1 train induces a second resonance on the bridge at approximately 222 km/h , close to the theoretical resonant velocity. The resonant velocity decreases with the soil stiffness, as a consequence of the reduction in the bridge fundamental frequency. In the case of the HSLM-A6 train, resonance does not take place because the span length and characteristic distance L∕d = 0.76 , is close to a cancellation condition. This phenomenon of cancellation of resonance occurs even in the case of the soil with c s = 150 m/s . From the figures it can be seen that when SSI is considered, the natural frequencies of the bridge computed following the present approach are very close to the results from the BE-FE model in Doménech et al. (2016). Nevertheless, the overall damping exhibited by the bridge-foundation-soil system is higher, especially for soils with very low values of c s .

Case study
In this section the dynamic response of Old Guadiana bridge (Fig. 12), from Madrid-Alcázar de San Juan-Jaén conventional railway line in Spain, is reproduced and compared with experimental measurements. It is a SS bridge with two identical spans of approximately 12 m of length. The deck is composed of two single track independent decks of 5.075 m width. Each deck consists of a 25 cm thickness concrete slab resting on five prestressed concrete beams with a 75 cm height rectangular cross-section (Fig. 13a). This bridge was monitored by the authors in 2019. A portable acquisition system LAN-XI of Brüel & Kjaer was used. The acquisition system fed the sensors (accelerometers) and an instrumented impact hammer in the case of the soil tests. LAN-XI also performed the Analog/Digital conversion (A/D). The A/D was carried out at a high sampling frequency, f s = 4096 Hz , that avoided aliasing effects. The acquisition equipment was connected to a laptop for data storage. Endevco model 86 piezoelectric accelerometers were used with a nominal sensitivity of 10 V/g and a lower frequency limit of approximately 0.1 Hz . Structural response recordings were decimated (order 16) to carry out data analysis in the frequency range of interest. Structural responses were filtered applying two third-order Chebyshev filters with high-pass and low-pass frequencies of 1 Hz and 30 Hz , respectively. The modal parameters of the bridges were identified from ambient vibration data by the stochastic subspace identification technique Fig. 12 a Old Guadiana bridge and b soil test set-up (Reynders 2012). The dynamic characterisation of the soil was carried out by the seismic refraction and the Spectral Analysis of Surface Waves (SASW) tests. More details of the experimental campaign, the identified properties of the bridge and surrounding soil and the dynamic response under several train passages may be consulted in Galvín et al. (2021).
The site has two upper soil layers with a shear wave velocity close to 300 m/s resting on a halfspace with c s = 250 m/s ( Fig. 13b and Table 1). According to Eurocode 8 (1998), the soil profile is approximated as a homogeneous stratum with average shear wave velocity V s30 = 303 m/s . Likewise, a material damping ratio of 7.85% is admitted. Following, the results obtained from the actual stratigraphy and the equivalent homogeneous soil are compared.
A first vertical bending mode (see Fig. 14) with a natural frequency of 9.8 Hz is identified from ambient vibrations. The modal damping corresponding to this mode obtained under traffic conditions reaches 2.8% , which is higher than the value prescribed by standards for design purposes for this particular length and bridge typology (1.5% as  Hz . This frequency reduces to 9.8 Hz due to SSI effects. Figure 15 shows the corresponding numerical mode shape which is the only bending mode of the beam bridge below 30 Hz . The MAC value (Allemang and Brown 1982) between the experimental mode shape corresponding to the longitudinal middle section and the numerical mode shape is 0.99. As part of the initial approach, only the contribution of this mode is taken into account as it dominates the vertical response at mid-span along the centre line. Nevertheless the first torsion and transverse bending modes with frequencies 11.0 Hz and 12.8 Hz were also identified. The response of the bridge is presented under the circulation of the Renfe Altaria train with 9 coaches crossing the bridge at v = 160 km/h . Table 2 summarizes the train coach distribution and axle distances. Figure 16 shows the bridge acceleration response at the mid-span section, measured under the track and its comparison with the numerical predictions obtained with the proposed methodology. The agreement between numerical and experimental results is good given the simplicity of the beam model used for the deck. The response contains peaks associated with the excitation (i.e. ratio of train speed v to the characteristic distance v∕d = 3.38 Hz ), corresponding harmonics and with the bridge lowest natural frequency. The experimental response shows certain participation of modes between 10 and 20 Hz that cannot be reproduced with the beam model. The train speed is close to a third resonance of the fundamental mode of the bridge: v 1,3 = f 1 × d∕3 = 9.8 × 13.14∕3 × 3.6 = 154 km/h . This justifies the dominant contribution of this mode in the frequency domain. The numerical prediction overestimates the response at the resonant frequency, which is expected because dissipation mechanisms such as vehicle-structure and track-structure interaction are not included, and are particularly relevant under these conditions. In this case, the actual stratigraphy is well represented by the equivalent homogeneous soil. The comparison between both results presents very small differences in the frequency range analysed.  The time used by a laptop computer with 1.9 GHz processor and 16 GB of RAM in order to solve the complete problem is 26.69 s : the homogeneous soil-foundation problem required 24.02 s (two frequencies were computed in this case), the eigenvalue problem required 0.92 s , and the modal superposition required 1.75 s , for one train speed. A speed sweep from [150∶1∶400] km/h was solved in 439.2 s for this particular train and soil-foundation-bridge system. This is considerably lower than the time required if the same problem is solved using BE-FE methods in the time domain.

Conclusions
A fast numerical methodology to determine the dynamic response of railway bridges under the passage of running trains including the effect of soil-structure interaction is proposed. The coupled problem is sub-structured into a soil-foundation and a soil-foundation-bridge interaction problems. In the first step, the foundation and surrounding soil are simulated via a FE-PML model. Using this approach, equivalent frequency dependent dynamic stiffness and damping characteristics of the soil-foundation system are computed. In the second step, the dynamic response of the structure under railway traffic is calculated using a FE model that includes spring/dashpot elements at the supports, using the equivalent properties extracted from the first problem. This model is solved by complex modal superposition considering the equivalent dynamic stiffness and damping of the soil-foundation corresponding to each natural frequency. The complex modal parameters are obtained iteratively and the dynamic complex modal equations are solved using a Runge-Kutta explicit algorithm. The benefits of the proposed approach increase dramatically with the problem complexity since it allows exploiting dedicated analysis tools for addressing the soil-foundation and the superstructure domain.
From the analyses performed the following is concluded: • The proposed approach can be used to obtain the dynamic response of railway bridges under running trains considering SSI, with minimal computational effort. However, it is remarked that soil-foundation impedances should be previously calculated. • The results obtained in terms of the bridge deck vertical displacements and accelerations are similar to those given by BE-FE analyses. • Resonance and cancellation phenomena, which are often responsible for the maximum bridge acceleration, are accurately reproduced. • In this particular study, the structural model is simplified to a beam bridge and superficial infinitely rigid plate foundations, because the main objective of the work is to present the method and to compare with previous results. However, the consideration of more complex geometries for the bridge and the foundations, and the simulation of stratified soils is also straightforward to implement. In particular, the effect of the foundation mass should be further analysed. • Using the proposed approach it is possible to solve the dynamic problem taking into account high modal contributions of the structure, with minimal computational cost. • The results are compared with those given by a fully coupled 3D BE-FE model solved in the time domain, where the Green's function for an elastic half-space is used as the fundamental solution for the soil, and material damping is assumed to be constant. In this case similar resonant speeds are determined, and close to the resonant frequency, a more highly damped response is obtained using the present method, particularly for soils with very low shear wave velocities. • Using the proposed model the experimental response of a real beam-type bridge under train loading is reproduced, capturing the frequency contributions of the first structural bending mode and the loading frequency contributions.