Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings

ANALYSIS AND MITIGATION OF VIBRATIONS INDUCED BY THE PASSAGE OF HIGH-SPEED TRAINS IN NEARBY BUILDINGS JOÃO MANUEL DE OLIVEIRA BARBOSA 2013 Supervisors: Prof. Álvaro Azevedo (FEUP) Prof. Rui Calçada (FEUP Prof. Eduardo Kausel (MIT) Thesis presented to the Faculty of Engineering of the University of Porto for the Doctor Degree in Civil Engineering Ao titio e à titia Aos meus pais Aos meus irmãos Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings i Abstract The present dissertation addresses the subject of vibrations induced by the passage of high-speed trains. The main objective of the study is the development of numerical tools that allow investigating distinct geometries of train and track, considering buildings in the proximity of the track, and also assessing the behavior of mitigation solutions. The problem of vibrations induced by moving vehicles is divided into three stages: a generation stage, in which the vehicle interacts with the track; a propagation stage, in which the forces that the train transmits to the track originate waves that propagate through track and ground; and a reception stage, in which the waves reach a nearby building, causing it to respond dynamically. Since the geometric specifications of the problem vary within the three stages, different strategies are chosen for each stage: 1. The generation stage involves a discrete structure (the vehicle) moving on top of a structure whose longitudinal dimension is infinite (track-ground system). In this way, the problem is formulated in a moving frame of reference, being the equations solved in the frequency domain; 2. For the propagation stage, since it is assumed that track and soil are invariant in the longitudinal direction, then the problem is formulated in the wavenumber-frequency domain (2.5D). In this way, the three-dimensional problem is reduced to a series of two-dimensional problems of smaller dimensions that are faster to solve. The track is simulated with finite elements while the surface of the soil interacting with the track is simulated with boundary elements. Mitigation measures in the soil must be included in this stage; 3. In the reception stage, the three-dimensional structure to be analyzed is irregular in all directions and therefore the 2.5D procedure cannot be applied. For this reason, a 3D frequency domain formulation is used, in which the structure is simulated with finite elements and the soil is simulated with boundary elements. The exterior loads considered in this stage are calculated based on the results of the propagation stage. The response of the soil, which is of great relevance for the problem, is accounted for through the boundary element method (BEM). The fundamental solutions used to nurture the BEM are obtained with the thin-layer method, being this the main difference between the strategy adopted herein and the procedures followed by other authors that also use the BEM. The numerical procedures mentioned above are described in chapters 2-4. Additionally, in chapter 4, the link between the three stages and the distinct procedures is established. In chapter 5, the methodology is applied to the study of trenches as mitigation solutions. Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings iii Resumo Esta dissertação debruça-se sobre o tema de vibrações induzidas por tráfego ferroviário de alta velocidade. O principal objetivo do estudo reside no desenvolvimento de ferramentas numéricas capazes de simular o fenómeno e que permitam investigar vias e comboios com geometrias variadas, considerar edifícios na proximidade da via, e ainda avaliar o comportamento de diversas medidas mitigadoras. O problema de vibrações induzidas por veículos pode ser dividido em 3 fases: uma fase de geração, em que o veículo interage com a via; uma fase de propagação, em que a ação que o veículo transmite à via-férrea origina ondas que se propagam através desta e através do solo; e uma fase de receção, em que as ondas chegam a um edifício próximo da via, induzindo a resposta dinâmica da estrutura. As especificidades geométricas variam de fase para fase, pelo que são tomadas diferentes estratégias para cada fase: 1. A fase de geração envolve a interação entre uma estrutura móvel de caráter discreto (o veículo) e uma estrutura de dimensão longitudinal infinita (sistema solo-via). Desta forma, formula-se o problema com base num referencial móvel, sendo as equações posteriormente resolvidas no domínio da frequência; 2. Para a fase de propagação, uma vez que se assume que tanto a via como o solo apresentam invariância longitudinal, as equações são formuladas no domínio do número de onda e da frequência (2.5D). Desta forma, reduz-se o problema tridimensional a um somatório de problemas bidimensionais de menor dimensão e de mais rápida resolução. A via-férrea é modelada com recurso a elementos finitos e a superfície do solo em contacto com a via é modelada com elementos de contorno. De salientar que medidas mitigadoras no solo, como por exemplo, trincheiras, devem ser incluídas nesta fase; 3. Na fase de receção, a estrutura a analisar apresenta um caracter tridimensional e irregular em todas as direções, pelo que a formulação 2.5D deixa de ser válida. Assim, opta-se por uma formulação tridimensional, no domínio da frequência, em que a estrutura é modelada com elementos finitos e a superfície do solo em contacto com a estrutura é modelada com elementos de contorno. Convém ainda referir que a ação considerada nesta fase é calculada com base nos resultados da fase de propagação. Relativamente ao solo, elemento preponderante em todo o problema, o seu comportamento é tido em conta por intermédio do método dos elementos de contorno (MEC). As soluções fundamentais usadas para alimentar o MEC são obtidas pelo thin-layer method, sendo esta a principal diferença entre a estratégia aqui adotada e a estratégia seguida por outros autores que também usam o MEC. A formulação das ferramentas acima mencionadas é descrita nos capítulos 2 a 4, sendo ainda no quarto capítulo exemplificada a ligação entre os diferentes procedimentos. No capítulo 5, a metodologia é aplicada ao estudo de medidas mitigadoras sob a forma de trincheiras. Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings v Acknowledgments It took me five years to conclude this Doctoral work, and such would not have been possible without the contribution and support of the people that one way or the other accompanied me throughout this process. To them, I want to demonstrate my appreciation. • To my supervisors, Professors Álvaro Azevedo, Rui Calçada and Eduardo Kausel, for having accepted the responsibility of guiding. I would like to stress the role of Professor Eduardo Kausel, who received me in the United States and who was essential for the success of my stay at MIT. • To Professor Pedro Costa, for his comments and opinions given during the development of the thesis, but mostly during the last year of work, which relates to Chapter 5 herein presented. • To the colleagues developers of FEMIX, Prof. Álvaro Azevedo and Sérgio Neves, for making the software available, and for their effort in making FEMIX an attractive choice for people who want or need to develop their own calculation routines. • To colleague Ricardo Nobre, PhD student in the field of computer science at FEUP, to the research group SPECS, to Professor Pedro Sobral from University Fernando Pessoa, to Professors João Cardoso, Rui Rodrigues and Jorge Barbosa from the Department of Computer Science of FEUP, and to Nuno Subtil from NVIDIA, for their support regarding GPU programming. • To all colleagues from the High-Speed research group, headed by Professor Rui Calçada, for all help provided during these five years. • To all colleagues and friends that have worked in office H301 and that made the working environment very friendly — Alés de Miguel, Cristina Ribeiro, Fernando Bastos, Luís Martins, Mário Marques, Miguel Araújo, Nuno Santos, Ricardo Monteiro, Sérgio Neves. • To the MIT office mates, for the same reason — Dilip Thk, Swapnil Rajiwade, Rositta Jünemann, Anthoula Agn, Zeid Alghareeb, Inez Azaiez. • To Fundação para a Ciência e Tecnologia (FCT – Portuguese Foundation for Science and Technology), that provided me financial support through PhD grant SFRH/BD/47724/2008 and through the research project PTDC/ECM/ 114505/2009 — Ground Vibration and Noise Induced by High-Speed Trains: Prediction and Mitigation. • To MIT-Portugal program, for facilitating my research periods at MIT. At last, the most special acknowledgment goes to those who have followed me on a daily basis, ever since my birth. An immeasurable thank you to my parents, Jorge Barbosa and Maria José Barbosa, my brothers, Ricardo Jorge and Rui Vasco, to may granduncle and grand aunt, Manuel Sala (titio) and Maria Elisa (titia). Without their absolute support, I would not have started the PhD. Chapter 1 – Introduction 2 number of complaints received from inhabitants who sensed the vibrations and feared damage on their properties. The construction of the subway systems inside the cities (excessively close to residences), together with the increase of train capacities and comfort standards, gave more importance to the problem of vibrations induced by moving vehicles and led to several experimental campaigns (Wilson et al., 1983; Dawn and Stanworth, 1979; Melke and Kramer, 1983). The experimental campaigns were carried out to assess if the vibrations induced by the circulation of trains could damage buildings, cause discomfort to inhabitants or cause the malfunction of sensitive equipment and, at the same time, to understand the mechanisms of generation of vibrations and to study possible mitigation measures. In more recent years, in part due to the continuous increase of the weight and speed of trains and in part due to the higher standards for comfort, the problem became even more important and, in addition to the problem of vibrations induced on buildings, a new problem emerged: at some lines resting on soft soils the train speed approached the propagating velocity of the waves and, as a consequence, large displacements in the embankment were observed, causing the risk of derailment of the train. This happened, for example, in the very well documented case of Ledsgard, Sweden (Hall, 2000), and in a line of the Northwest of France (Picoux and Le Houedec, 2005). In those lines, the circulation speed was reduced and, in several other lines around the world, new experimental campaigns were conducted to assess the vibration levels. In addition, studies were performed on countermeasures to mitigate the vibrations and prediction models were developed or improved. Two strategies for the study of the phenomenon of vehicle induced vibrations can be adopted: field measurements and numerical/analytical predictions. Field measurements are performed under real conditions and provide an enriched set of results, which account for all the factors that influence the phenomenon. The results of the experimental campaigns indicate which aspects most influence the level of induced vibrations, thus showing which factors must be taken into account when developing numerical models. When large databases collecting experimental data are available, it is possible to perform extrapolations in order to predict the vibration levels for scenarios with similar conditions (soil, track, vehicle, structure), and, consequently, it is possible to develop empirical models. The drawback of field experiments is their high cost and so it is desirable that they are employed as less as possible. Nevertheless, experiments are always needed to obtain inputs for numerical models and also to validate these models. On the other hand, to use prediction models is not as expensive as performing in situ experiments but, contrarily to experiments, these models cannot reproduce the whole reality since they are based on assumptions and simplifications. On one side, empirical models that are developed based on experimental data provide good estimates but their range of applicability is limited to scenarios with conditions similar to the experiments. On the other side, numerical models are versatile and allow studying the influence of certain parameters on the vibration levels, but their accuracy depends on the simplifications made and on the assumptions on which the models lie. Numerical models can usually be adapted in order to account for countermeasures and study their performance. In the following sub-sections, a historical overview of experimental campaigns, prediction models and countermeasures for vibrations is presented. Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 3 1.2.1 Experimental campaigns As mentioned earlier, the problem of vibrations induced by moving vehicles harks back to the XIX century, when South (1863) performed experimental investigations in a tunnel located in Watford, England. These experiments were carried out to sustain the idea that the passage of freight trains in Greenwich Park, near the Royal Observatory, could be harmful to its equipment. The Watford tunnel was selected for the experiments because the geotechnical conditions at that place were similar to the ones observed in Greenwich. Vibrations induced by road traffic have also been subjected to investigations since the mid twentieth century, as reported by Sutherland (1950), who refers to experiments made in Canada. There, the inhabitants of the city of Winnipeg sensed the vibrations and feared damage on their properties, and therefore presented their complaints to the City Hall. As a consequence of the multiple complaints, the National Research Council of Canada promoted experiments to investigate the problem. The results of the experiments suggested that the induced vibrations, even if felt, were not strong enough to cause damage in the structures. In addition, it was concluded that the irregularities of the surface of the road (such as bumps) were the factor that contributed most to the vibration levels. Also concerning road traffic induced vibrations and still in Canada, more experimental campaigns were carried out to understand the phenomenon. As a consequence of complaints lodged by residents in Quebec, Al-Hunaidi and Rainer (1991a) performed field experiments in two different sites in order to study the factors that influence the level of vibrations. They studied the influence of speed, weight and type of vehicle and road roughness, and concluded that while the vehicle speed and the road irregularities influenced considerably the vibrations, the mass of the vehicle had little influence. Later, Al-Hunaidi and collaborators performed more tests in nine different sites of the city of Montreal to study the influence of the suspension system of the vehicles (Al-Hunaidi et al., 1996) and concluded that the vibrations could be greatly reduced by imposing an axle hop frequency lower than the cutoff frequency of the soil. The influence of the road surface condition and seasonal variation of soil conditions has also been taken into account in another study (Al-Hunaidi and Tremblay, 1997). To study the efficiency of soil improvement as a countermeasure for road traffic induced vibrations, Taniguchi and Okada (1981) measured the acceleration before and after the improvement of the soil using the lime pile technique. When they compared the spectrum of the vibration reduction with the acceleration spectrum of a point situated 8m away from a national road, they observed that the frequency range over which the vibrations were reduced due to the soil improvement was approximately the same range that was excited by the road traffic, thus concluding that the method was efficient. To validate their numerical model for the prediction of road traffic induced vibrations, Lombaert and Degrande organized two campaigns where a truck moving with variable speed was submitted to an artificial unevenness (Lombaert and Degrande, 2001; Lombaert and Degrande, 2003). Free field vibrations and acceleration of the axles were measured and the dynamic properties of the road and soil were determined experimentally. The results suggested that the vibration levels were dependent on the vehicle speed, the shape of the unevenness, and the vehicle and soil characteristics. The comparison between experimental and predicted results showed some differences that were explained by the loss of contact between the rear axle and the road. Chapter 1 – Introduction 4 With respect to rail traffic induced vibrations, Wilson et al. (1983) reported experiments made during the seventies to study the use of floating slabs as a countermeasure to mitigate groundborne vibrations and to study the influence of the properties of the bogies on the induced vibrations. Measurements were made in tunnels, on the free surface, and inside buildings. Later, in the late seventies, Dawn and Stanworth (1979) measured the vibration levels on a wall of a single storey building situated about 42m away from a track during the passage of trains at speeds up to 100 km/h. They noticed that the vibrations increased with the speed of the train and observed that the frequency content of the response presented a peak on the passage frequency of the sleepers. Dawn (1983) performed further experimental studies and confirmed that the passage frequency of the sleepers is indeed a mechanism of excitation. He also recognized that the critical velocity (to which corresponds the maximum ground response) occurred when the sleeper passage frequency coincided with the resonance frequency of the vehicle-track system. Melke and Kramer (1983) reached the same conclusions in their experimental studies. In more recent years, with the increase of the train speed, the problem of induced vibrations has been given even more importance and new experimental studies have been performed. In addition, with the construction of new railway lines, their homologation tests enabled the execution of new field experiments. In Germany, Auersch (1994; 2005) performed measurements at three different sites near Würzburg during test runs of the ICE train with different configurations and at speeds between 100 and 300 km/h. During the tests, which considered three different track conditions (surface line, bridge and tunnel), the vibrations of the vehicle, track and soil were recorded. The results showed that the quasi-static component of the axle load was important for the response of the track and the surrounding soil, and that its importance vanished rapidly with the distance. The results also suggested that the sleepers act as harmonic forces, whose intensity increases with the train speed, but remains constant when the sleeper passage frequency exceeds the vehicle-track resonance frequency. In Sweden, as a consequence of the high vibration levels observed shortly after the opening of the line between Göteborg and Malmö, in 1997, the train speed was reduced at some locations and investigations were conducted during the Autumn of 1997 and Spring of 1998 to diagnose the problem and to find solutions (Madshus and Kaynia, 2000; Hall, 2000). A X-2000 passenger train was used in a total of 20 runs at speeds ranging from 10 to 202 km/h and the responses of rail, sleepers, embankment and ground (at the surface and its interior) were measured. It was observed that for speeds below 70 km/h the displacements of the ground were similar to those obtained considering static loading and therefore were independent from the speed. At speeds around 200 km/h, the amplitude of the displacements increased drastically, causing the risk of derailment of the train. In Belgium, the expansion of the railway network allowed to experimentally investigate the phenomenon in newly built high-speed lines. In December of 1997, six weeks before the inauguration of the high-speed line between Brussels and Paris, an extensive experimental campaign was organized by the Belgian railway company during the homologation phase of the line. Track response and free field vibrations up to 72 meters away from the track were measured during the passage of a Thalys train at speeds varying between 223 and 314 km/h (Degrande and Schillemans, 2001). Five years later, in August and September of 2002, the high-speed line between Brussels and Köln was also submitted to homologation tests (Kogut et al., 2003). The tests were performed at two different sites, Lincent and Waremme, and Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 5 included the measurement of vibration levels at the track, in the free field and in a single family dwelling located 50m away from the track during the passage of Thalys trains and IC trains at variable speed. The tests were complemented with the experimental determination of some dynamic properties of the track and soil. The two sets of experiments showed that the passage frequencies of bogies and sleepers and their higher harmonics could be noticed on the spectrum of the responses, namely in the near field. Also, differences in the registered responses for the different trains suggested that the induced vibrations depend on the train properties. In both cases, attenuation of vibrations with the distance to the track was detected. In Italy, Lai et al. (2005) measured the transfer functions in two sections of a tunnel in the city of Rome. The aim of the study was to assess if the level of vibrations would affect the surrounding buildings. Since the line was not yet operational, it was not possible to perform a direct measurement of vibrations induced by the railway traffic. Hence, the transfer functions from the tunnel to the free field and to the interior of buildings were determined experimentally using a mechanical hammer. The transfer functions would serve as inputs in a simple numerical model to predict the level of vibrations induced by future passing trains. In the Northwest of France, after it was observed that the ground presented excessive displacements, Picoux and Le Houedec (2005) measured the vibrations on rails, sleepers and free-field during the passage of different trains. The results showed the influence of the train speed and of the type of train on the induced vibrations. In England, within the framework of the CONVURT project, vibrations were measured at a site in Regent’s Park, London, during 35 passages of a test train in a tunnel at speeds between 20 and 50 km/h (Degrande et al., 2006). Accelerations of the axle boxes, of the tunnel, of the free field (both at surface and inside the soil) and on several floors of two buildings situated 70m away from the tunnel were measured. Rail and wheel roughness have also been measured and track characteristics were determined by receptance tests. Analysis of the measured fields allowed concluding that the peak velocities on the axle boxes and track increased with the train speed, a tendency that was less pronounced in the free-field and in the buildings. In Beijing, China, a subway line was planned to pass close to the Physics Laboratory of the Beijing University and so there was concern about the vibrations that would be induced by the rail traffic (Gupta et al., 2008). To study if certain equipments would need to change place, measurements were performed in the free-field near the lab and inside the building to evaluate the existing vibration levels (induced by road traffic and people). In addition, measurements were made in a different line of the Beijing subway system with similar characteristics. The superposition of the existing background and the predicted vibration levels would provide the vibration level expected in the labs. In Northeast China, Xia analyzed experimentally the problem of vibrations in buildings induced by trains running on bridges. Trains running at speeds varying between 60 and 80 km/h (Xia et al., 2005a) and between 160 and 307 km/h (Xia et al., 2005b) were considered. It was observed that the vibration levels would increase with the weight and speed of the train and would attenuate with the distance to the railway line. It was also observed that the vibrations were stronger at higher floors, exceeding in some places the levels allowed by the Chinese code. To finalize the field measurements, some experimental campaigns have also been organized at the Portuguese Northern line by researchers from FEUP, whose objective was to characterize the site conditions and measure the vibrations induced by real traffic (Alves Costa et al., Chapter 1 – Introduction 6 2012a; dos Santos, 2013). At the same site, further investigations have been conducted in order to evaluate the dispersion of the responses along the longitudinal directions. The results obtained from these campaigns are expected to be published in the near future. 1.2.2 Prediction models Two types of prediction models can be considered: empirical and analytical/numerical. Empirical models are based on results from experimental campaigns and usually provide good predictions, but their range of applicability is limited to scenarios that are similar to the conditions under which the experiments were performed, thus lacking versatility. Description of empirical models can be found in works by Kurzweil ( 1979), Melke (1988) and Madshus et al. (1996). These models use chains of transmission losses for the source-path-receiver system and consider parameters such as train speed, axle loads, suspension systems, weight of the train, wheels and rail conditions, rail fastening systems, type of track, type of tunnel and type of buildings. The model described by Madshus et al. (1996) was developed based on a large number of vibration measurements made in Norway and Sweden and was used for the planning of a high-speed railway line in Norway. In another empirical work, to evaluate problems of excessive vibrations in preliminary stages, Bahrekazemi (2004) presented a model that is based on measurements performed in several sites of Sweden. On the other hand, numerical and analytical methods are more versatile and can be efficiently used to study the effect of train speed or weight, track type, material resiliency, ground conditions, etc. The drawback is that these methods rely on idealizations and simplifications, then failing to reproduce reality as accurately as it would be possible with field experiments. Nonetheless, depending on the degree of detail of the model, the obtained prediction can be acceptable and useful. To correctly model vibrations induced by vehicles, three stages must be accounted for in a numerical/analytical model: the generation stage, the propagation stage and the reception stage (Figure 1.1). In the generation stage, the vehicle interacts with the track and induces a moving stress field on it. The stresses are transmitted from the vehicle to the track through contact surfaces (wheels or tyres) that move in space. Due to the dynamic behavior of the vehicle and its interaction with the track, the vehicle is subjected to accelerations and so the contact stresses, besides moving with the vehicle, also change their value with time. The non varying component of the contact stresses is called quasi-static excitation (forces per wheel or tyre) while the component varying with time is termed dynamic excitation. In the propagation stage, the stress fields (or the vibrations) propagate through the track and part of them is transmitted to the soil. These stresses continue to propagate in the soil, being reflected or refracted whenever a different material or a barrier is encountered, and finally reach the building. In the reception stage, the vibrations that reach the building induce a dynamic response on it. The problem of vibrations induced by moving vehicles is three dimensional: the vehicle moves in one direction while the waves propagate in the soil in three directions. Modeling a three dimensional problem can become very complicated and time consuming, even for the current computers. For this reason, the first works assumed that the phenomenon could be described by 2D models. For example, in their review paper, Gutowski and Dym (1976) mentioned that the vibrations generated along a road or a railway track could be modeled as a line source as long as the roadway was relatively uniform and the receiver was in the far field, but close enough to the source (less than 1/π times the length of the roadway or the train). The authors supplemented that if the ground motions were dominated by the surface (Rayleigh) Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 7 waves, then there would be no geometric damping and so the vibrations would attenuate only due to material damping. Later, Verhas (1979) compared the results of line source models with the results of point source models and concluded that by neglecting the geometric damping of waves inaccurate predictions would be obtained. This author suggested that the combination of the two models would yield better predictions, but no guidelines on how to combine the results from each model were indicated. In another 2D work where a finite element (FE) model was used, in order to account for the geometric damping of surface waves, the accelerations were corrected by a factor 1/ r , being r the distance to the source (Taniguchi and Okada, 1981). This methodology was used to study the efficiency of soil improvement via the lime pile technique as a countermeasure. Also using a 2D FE procedure, Chua et al. (1995) determined the vibration levels in a four-storey podium block due to the passage of trains in a double-box tunnel, accounting both for the quasi-static and for the dynamic excitation. The authors used an iterative nodal condensation procedure to avoid extremely large meshes. Figure 1.1: Generation, propagation and reception of vibrations (Hall, 2003) With the improvement of computational performance, both in terms of memory and speed, the use of 3D models became possible. One of the first works that considered the threedimensionality of the problem was performed by Krylov (Krylov and Ferguson, 1994; Krylov, 1994; Krylov, 1995). In their work, Krylov and his collaborators developed a model for surface trains where the forces transmitted to the ground through each sleeper are calculated analytically, and then, considering the sleepers as point sources, the field induced by each sleeper is combined, thus obtaining the response of the soil due to the passage of the train. Only the quasi-static component of the excitation is considered. The method for the calculation of the forces transmitted to the soil has been used by other researchers. Chapter 1 – Introduction 8 Conventional methods used for the analysis of three-dimensional problems, such as the FE method and the boundary element (BE) method, have also been employed to analyze the problem of induced vibrations. The FE method requires the discretization of the domain, which for 3D problems results in a large number of degrees of freedom and in sparse symmetric matrices. This method can be used to model irregular domains and, when applied in the time domain, can account for the non-linear behavior of materials. By itself, the classical FE method cannot simulate infinite domains, so special procedures need to be considered at the boundaries of the truncated domains in order to avoid fictitious reflections. Contrarily, the BE approach only requires the discretization of the boundary of the domain, thus resulting in less degrees of freedom, but, unlike the FE method, leads to full nonsymmetric systems of equations. The BE method takes into account the radiation of waves towards infinity, but cannot account for non-linearities and requires the knowledge of the so called Green’s functions (GF) or fundamental solutions. The hybrid FE-BE method combines the advantages of both approaches, being its use very attractive when the coupling between irregular domains and unbounded domains is required. Regarding the FE approach, Hall (2003) used a time domain methodology and treated the reflections at the boundaries using dashpots. The considered mesh led to reasonable results only up to the frequency 10Hz and in the calculations only the quasi-static component of the excitation was considered. The results of the model showed a transient phenomenon that was not observed in real measurements and that was originated by the entrance of the loads in the model. However, this numerical phenomenon would have dissipated due to damping by the time that the waves reached the other extreme of the model, and so the results at that extreme were better. Using a model with 65 meters in the longitudinal direction, good results were limited to the near field. To obtain better results farther from the track, longer models would be needed, which would render the mesh impractical for calculation. The same author compared the results obtained with 3D models with those obtained with simpler 2D models and concluded that the 2D models could be used to study certain effects of traffic induced vibrations but not to obtain good predictions of the induced levels of vibrations (Hall, 2000). In another work using the FE approach, Ekevid and Wiberg (2002) followed a similar approach but instead of treating the boundaries of the mesh with dashpots, these authors used the scaled boundary finite element method (Wolf, 2003). Even though the proposed methodology accounted for the radiation of waves to infinity, the fact of using a 3D mesh required a very large computational effort. Also following the FE approach, Ju used a 3D formulation to simulate soil vibrations due to a high-speed train crossing a bridge and to study the efficiency of trenches (Ju, 2002) and of soil improvement (Ju, 2004) as countermeasures. The boundaries were treated with first-order absorbing boundaries and the systems of equations were solved using the preconditioned conjugate gradients method, i.e., an iterative method. The calculation time of the problem was over one week. As for the BE and hybrid FE-BE approaches, Bode et al. (2002) used the BE method to model the soil and the FE method to model the sleepers and the rails (dos Santos, 2013). The methodology was formulated in the time domain and was used to determine the vibrations in the free-field and to study the influence of the soil-sleeper coupling scheme. The GF considered for the soil were the half-space Green’s functions, thus limiting its discretization to the regions interacting with the sleepers. A similar strategy was followed by O'Brien and Rizos (2005), but instead of using half-space GF, they used full-space GF, which demanded Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 9 the additional discretization of the free surface of the soil and, consequently, increased the number of degrees of freedom. In his PhD thesis, Galvín used an iterative scheme to couple the FE and BE methods and evaluated the response of structures near the railway track (Galvín, 2007). To account for the train excitation, he considered the forces transmitted by the sleepers to the soil or ballast as given by Krylov, and so only the quasi-static component of the excitation is included. Auersch (1994) presented a model for surface trains based on transfer functions of point loads that could be either determined experimentally or calculated numerically. The train was simulated by a chain of point loads representing the axles, and for each axle a force function was assumed with the intent of simulating the dynamic forces of the train. Later, the work was extended and the FE method was combined with the BE method in order to include the traintrack-soil interaction and to consider the irregularities of the vehicle/track system and the discrete sleeper support (Auersch, 2005). In more recent years, approaches that take advantage of the invariance or periodicity of the geometry in the longitudinal direction were developed. Some authors consider that the geometry is invariant in the longitudinal direction, and after performing a Fourier transform of the field of variables in that direction, they reduce the three-dimensional problem to a series of two-dimensional problems. The transformed problems are solved in the wavenumberfrequency domain, which is termed 2.5D domain. Other authors consider that the problem is periodic in the longitudinal direction and after performing a Floquet transform of the field of variables in that direction, they reduce the geometry of the problem to a reference cell. Even though these two approaches, by themselves, cannot be used to predict the vibrations inside buildings (buildings are not invariant nor repeat themselves till infinite), they can be used to determine the wave fields that reach the buildings. Those wave fields can later be employed in the calculation of the response of the buildings. This is the approach followed by François, who obtains the response of buildings submitted to an incoming wave field using a FE-BE method in the time domain and considering the non-linear behavior of materials (François et al., 2006; François, 2008), and by Fiala et al. (2007), who use instead a frequency domain approach. The following group of works assumed invariant geometries and used the 2.5D approach. Dieterman and Metrikine (1996) coupled a beam and a half-space to model the track-soil interaction problem. They assumed smooth contact between the beam and the soil (no transmission of shear stresses) and assumed a uniform distribution of normal stresses between the beam and the half-space over the width of the beam. The model was used to determine the critical velocities of moving loads on the track-soil system. The authors concluded that there were two critical velocities: one that corresponds to the Rayleigh wave velocity and the other being slightly smaller. Both velocities resulted in severe amplifications of the beam displacements. Metrikine and Popp (2000) solved the same problem considering a viscoelastic layer instead of a half-space. They concluded that the critical speed of the moving load was close to the Rayleigh velocity of the layer and that it decreased slightly with the depth of the layer. The critical velocity of harmonic loads was treated by Dieterman and Metrikine (1997). These authors simulated the ballast using an elastic layer and determined the speed of the harmonic load that caused resonance of the system as a function of the thickness of the layer and of the frequency of the load. They concluded that resonance occurred when the velocity of the load was equal to the group velocity of the waves generated by the load. In the work by Steenbergen and Metrikine (2007), the validity of the assumptions concerning the Chapter 1 – Introduction 10 contact between the beam and the soil is studied and the authors concluded that in general, as long as the wavelength is large compared with the width of the beam and as long as the moving load spectrum does not present high frequency components, the simplified assumptions succeed in obtaining the track response. However, in order to obtain the near field response to constant moving loads and the far field response to moving harmonic loads with high oscillating frequencies, the interface between the soil and the track must be adequately modeled. In a similar line of investigation, Jones, Sheng and Petyt modeled the track as a layered beam (accounting for rail, railpads, sleepers and ballast) and the ground as a layered half-space, for which they used the transfer matrices derived by Haskell (1953) and Thomson (1950). The contact between the layered beam and the layered half-space follows the same assumptions as the work of Metrikine and collaborators. The model can account for fixed harmonic loads (Sheng et al., 1999a) and moving loads with constant or oscillating amplitude (Sheng et al., 1999b; Jones et al., 2000). Together with the forces transmitted by the train to the track (Jones and Block, 1996; Sheng et al., 2004) these models can be efficiently used to simulate the vibrations induced by the circulation of surface trains at variable speed, accounting both for the quasi-static and the dynamic components of the excitation. In the work by Sheng et al. (2003) this model is validated against the results of experimental measurements, and later it is used to study the influence of the track stiffness and the layered soil properties (Sheng et al., 2004). Still considering a beam resting on a layered half-space, Lombaert et al. (2000) developed a numerical model for road traffic induced vibrations. The road is modeled with a beam, the ground is modeled using the boundary element method, and the contact between the two substructures is assumed to be smooth. Since the tyres are much more flexible than the track-soil system, the vehicle-track interaction is uncoupled from the rest of the problem and the dynamic component of the excitation is calculated by simply submitting the vehicle to an irregular profile. The methodology is validated by means of field tests in the works of Lombaert and Degrande (2001, 2003) and used by Lombaert et al. (2001) to study the influence of the soil stratification. Clouteau et al. (2001) extended the methodology to the case of rail traffic and used it to exemplify the dynamic behavior of concrete slab tracks, to study resilient materials under the rail and slab, and to study the influence of the soil stratification. This last model was validated against experimental results (Lombaert et al., 2006a) and used to study the behavior of floating slabs as control measure for ground borne vibrations (Lombaert et al., 2006b). The model was also used to study the influence of the quasi-static and dynamic components of the excitation (Lombaert and Degrande, 2009). Also with respect to surface trains, Karlstrom and Bostrom (2006) developed a semianalytical model in which the ground is modeled as a visco-elastic layered half-space, the embankment and ballast with a rectangular layer, the sleepers by means of an anisotropic Kirchhoff plate, and the rails with Euler-Bernoulli beams. A simplification is considered at the lateral surfaces of the embankment, where it is assumed that the tangential stresses and the normal displacements are null. The results obtained with the model proposed by the authors were compared with the results obtained with finite element models, and it was concluded that the simplification provided good results as long as the load was vertical. The model was used to study the effect of the acceleration and deceleration of trains (Karlstrom, 2006), and it was concluded that the differences in terms of the vertical displacements between a train moving at constant speed and a train accelerating or decelerating were very small. In the longitudinal direction, however, large differences could be observed. These works were further extended Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 11 in order to account for the presence of water in the ground (poroelasticity) and to study the effect of trenches in the isolation of vibrations (Cao et al., 2012). As final example concerning surface trains, Alves Costa (2011) used 2.5D finite/infinite elements to simulate the Leedsgard case taking into account the large deformations that may exist due to the poor quality of the soil. With that purpose, he simulated the non-linear behavior of the soil under the track by considering equivalent elastic parameters that depended on the deformation level of the finite elements (Alves Costa et al., 2010). Additionally, he also used a 2.5D coupled FE-BE model to simulate the vibrations induced by the passage of trains in the Portuguese railway system (Alves Costa et al., 2012a). This last model was also used to study the strategy for modeling the train (Alves Costa et al., 2012b), and to study ballast mats as mitigation measures (Alves Costa et al., 2012c). From these two works it is concluded that train models can be reduced to axles, bogies and primary suspension systems, and that ballast mats perform better when placed beneath the subballast layer, and not as well when placed between ballast and subballast layers. Regarding tunnels, Forrest and Hunt (2006b) developed a model that assumes a tunnel with a cylindrical shape surrounded by a soil of infinite extent. Since the soil is treated by means of wave equations of an elastic continuum, no free surface is considered and consequently no surface waves are excited. In the far field, it is likely that buildings receive more energy from such waves than from body waves. Nonetheless, the model can be very effective in the evaluation of the response near the tunnel, where surface waves have much less influence. The model was also extended to account for tracks and then was used to assess the behavior of floating slab tracks (Forrest and Hunt, 2006a). The authors concluded that floating slabs yielded modest insertion losses and that under certain conditions they could even increase the transmission of vibrations. In the work by Hussein and Hunt (2007), the model was further extended to account for tangential forces at the tunnel walls, and in the work by Hussein et al. (2008) it was improved in order to permit the analysis of tunnels embedded in layered halfspaces. For that purpose, the authors assume that at a first step, the near field displacements are controlled by the dynamics of the tunnel and of the surrounding layer, i.e., they neglect the contribution of the other layers. At a second step, the response of the far field is calculated using the tractions calculated during the first step and assuming the proper stratification of the soil. Yang et al. (2003) approached the problem using finite/infinite elements formulated in the wavenumber-frequency domain (Yang and Hung, 2001). When compared with the 2.5D models mentioned so far, this approach has the advantage of considering the transverse stiffness of the track and of allowing more complex geometries both for track and soil. However, it may become less attractive because the number of dofs increases significantly. The authors used the developed methodology to study the stiffness, damping and stratification of the underlying soil, and concluded that increasing the stiffness results in the decrease of vibration levels, that increasing the damping results in a decrease of vibration levels only if the loads move faster than the Rayleigh wave velocity of the soil, and that the soil stratification is extremely relevant owing to the fact that the cut-off frequencies depend on the layers depths and because no waves can propagate below the first cut-off frequency. Also regarding tunnels, Rieckh et al. (2012) developed an invariant model in which the anisotropy of the soil is considered. The model is based on a boundary element formulation, and uses the 2.5D fundamental solutions of layered and anisotropic media, which are calculated with the method of potentials. Chapter 2 – Wave propagation in the soil: fundamental solutions 18 The four works referred to above were developed based on the Cargniard-de Hoop technique (De Hoop, 1960), which allows the direct evaluation of the double integral needed to transform the displacements from the wavenumber-frequency domain to the space-time domain. Interestingly, the direct evaluation of only one of these integrals is not possible and so closed form expressions of the fundamental solutions of half-spaces exist only in the spacetime domain. Such is so because when using contour integration to evaluate the improper integrals to transform the displacements from the wavenumber to the space domain (Erigen and Suhubi, 1975) or from the frequency to the time domain (Park and Kausel, 2004a), some branch integrals are obtained and these cannot be solved analytically. For the case of layered domains or sources/receivers inside homogeneous half-spaces, no closed form expressions are available and therefore it is needed to resort to numerical tools to determine the corresponding fundamental solutions. The most commonly used tools are based on integral transformation techniques, in which the fields of displacements are transformed to the wavenumber-frequency domain and consequently the wave equations are solved in that transformed domain. When necessary, the displacements can subsequently be transformed back to the space domain and/or time domain through the numerical evaluation of the integrals that result from the inverse transformations. In the transformed domain, the solutions can be found using the transfer matrices derived by Thomson (1950) and corrected by Haskell (1953), using the stiffness matrices derived by Kausel and Roesset (1981), using the method of Potentials (Tadeu et al., 2001; Tadeu and Antonio, 2001; Tadeu and António, 2002) or using the Thin-Layer Method (TLM) (Kausel and Peek, 1982). The first three approaches handle the propagation of waves within each layer without any approximation. In opposition, the TLM approach is based on discretizations of the domain in the vertical direction and in approximations of displacements within the layers by means of the interpolation functions. Its advantage over the other three approaches is that it enables the analytical evaluation of at least one inverse transformation. As a drawback, it requires the solution of two eigenvalue problems. In this work, the TLM is adopted for the calculation of the fundamental solutions of layered soils, being the procedure and its 2.5D formulation presented in this chapter. 2.2 Thin-Layer Method in Cartesian coordinates The TLM was introduced in the seventies (Lysmer, 1970; Lysmer and Waas, 1972; Waas, 1972) and since then it has found use in several areas related to wave propagation in layered media and in soil-structure interaction problems. The TLM is a semi-discrete numerical technique used for the analysis of wave motion in layered media, and consists in a finite element discretization in the direction of layering (for the case of soil, the vertical direction) combined with analytical solutions for the remaining directions, along which the material properties are assumed to be constant. A brief historical description of the method can be found in Park (2002). With respect to wave fields induced by moving loads or vehicles, the TLM is used in the works by Hanazato et al. (1991), Jones and Hunt (2011, 2012) and Celebi and Schmid (2005). In the first three references, the TLM is used in the context of transmitting boundaries (Kausel, 1988) or super-elements (Tassoulas and Kausel, 1981). For the presentation of the TLM, the wave equation is first expressed in matrix notation and is then discretized in the vertical direction. These steps follow the works from Park (2002) and Barbosa and Kausel (2012). Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 19 One note before starting: cross-anisotropic materials are more general than isotropic materials since they can reproduce different behaviors of the medium when loaded in different directions. For the case of the soil response, this is an important aspect because the layers are formed by vertical sedimentation of particles and therefore they present different mechanical characteristics in the vertical and horizontal directions. The drawback of considering crossanisotropic materials is that they require the quantification of more elastic constants, thus requiring more experiments in order to obtain the corresponding parameters. The constitutive matrix D of a cross-anisotropic material is the positive definite matrix defined by ( ) 2 2 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 t t t t t t t t t t G G G G D G G G D G G λ λ λ λ λ λ λ λ λ λ λ + >     + >     + > =  + >           D (2.1) where λ and G are the Lamé constants in the isotropic plane (horizontal planes) and t λ , t G , t D are the Lamé constants and the constrained modulus in the transverse direction (vertical direction). When t λ λ = , t G G = and 2 t D G λ = + , the material reduces to an isotropic one. With the intention of being more general, the following formulation considers crossanisotropic materials. Consider a horizontally homogeneous and vertically stratified cross-anisotropic elastic medium of infinite lateral extent and characterized by the depth dependent mass density ρ and the depth dependent constitutive matrix { } ( ) , 1,...,6 ij d i j= =D as defined in equation (2.1). Assume that the medium is subjected to an arbitrary dynamic load b placed at some location. With dots denoting partial derivatives with respect to time, the dynamic equilibrium equation at any point can be written compactly in matrix format as T ρ − = u L σb ɺɺ (2.2) where the displacement vector u , the stress vector σ and the differential operator L are defined as T x y z u u u   =   u (2.3) T xx yy zz yz xz xy σ σ σ σ σ σ   =   σ (2.4) T 0 0 0 0 0 0 0 0 0 x z y y z x z y x   ∂ ∂ ∂   ∂ ∂ ∂     ∂ ∂ ∂ =   ∂ ∂ ∂     ∂ ∂ ∂   ∂ ∂ ∂   L (2.5) Additionally, consider the stress-strain and strain-displacement relations = σ D ε (2.6) Chapter 2 – Wave propagation in the soil: fundamental solutions 20 = ε Lu (2.7) T xx yy zz yz xz xy ε ε ε ε ε ε   =   ε (2.8) The substitution of equations (2.6) and (2.7) in equation (2.2) results in the elastic wave equation (in the 3-D space) T ρ − = u L DLu b ɺɺ (2.9) The differential operator L can be expressed as xyz x y z ∂ ∂ ∂ = + + ∂ ∂ ∂ L L L L (2.10) where the matrices x L , y L and z L are 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 000 001 010 0 0 1 0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 x y z                         = = =                                     L L L (2.11) Since the domain under study consists of homogeneous horizontal layers, the material properties are piecewise constant with depth and invariant in the horizontal directions leading to ( ) 0 , , x y z α α ∂ ∂ = = D . Thus, the term T L DL in equation (2.9) can be expanded to ( ) ( ) ( ) 2 2 2 T 2 2 2 2 2 2 xx xy yx xz zx yy yz zy zz x x y x z y y z z ∂ ∂ ∂ = + + + + + ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ + + + + ∂ ∂ ∂ ∂ L DL D D D D D D D D D (2.12) being the material matrices αβ D defined by T , , , , x y z αβ α β α β = =D L D L (2.13) and given in Appendix I. Now, consider the internal stresses in horizontal planes, which are calculated by TT T zx zy zz z z σ σ σ   = = =   s L σL DLu (2.14) If one removes any horizontal slice of the medium and treats it as a free body in space, the dynamic equilibrium dictates the need to balance the internal stresses at the now exposed upper and lower surfaces with the external tractions t , i.e., u u l l     = =     −     t s t t s (2.15) where u t and l t are the external tractions applied at the upper and lower boundaries of the removed domain and u s and l s are the internal stresses at the same locations. Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 21 The first step for the formulation of the TLM is to discretize the domain in the vertical direction, i.e., to subdivide the medium into horizontal layers which are thin in the finite element sense, or in other words, which are small in comparison with the expected wavelengths and strain gradients. Thereafter, considering an arbitrary thin-layer as a free body in space (Figure 2.1), the displacements field inside the layer is approximated by means of interpolation functions, i.e. = u NU (2.16) where ( ) , x y =U U is a vector containing the nodal displacements (the nodes represent horizontal surfaces) T T T T 1 ... , , 1,2, , m j xj yj zj u u u j m     = = =    U u u u ⋯ (2.17) and ( ) z =N N is an interpolation matrix of the form [ ] 1 ... m N N= N I I (2.18) with j N being the interpolation functions, which depend on the vertical coordinate z , and with I being a 3 3 × identity matrix. (The subscript m is the number of nodal surfaces in each thin-layer, and 1 m − is the interpolation order. When 2 m > , there exist inner surfaces that are equidistant from each other. For example, 3 m = corresponds to a quadratic interpolation with one internal nodal surface, as shown in Figure 2.1, which depicts one thin-layer as a free body in space, acted upon and dynamically equilibrated by appropriate tractions applied onto the nodal surfaces.) Figure 2.1: Discretization into thin-layers and thin-layer as a free body in space ( 3 m = ). When substituting the interpolation (2.16) into the wave equation (2.9) and boundary conditions (2.15), it can be verified that these equations are not satisfied exactly because the interpolation is only an approximation of the actual field. As a result, one finds unbalanced body forces r and boundary tractions q of the form T ρ − + = b u L DLu r ɺɺ (2.19) 1 1 1 m m m       − = =       −       t s q q t s q (2.20) The discrete wave equation is obtained by applying the method of the weighted residuals and by requiring the virtual work done by the unbalanced forces within the thin-layer and on its bounding surfaces to be zero. This results in the discrete thin-layer equation h 1 1 , u t , m m u t x x y y z z Chapter 2 – Wave propagation in the soil: fundamental solutions 22 2 2 2 2 2 xx xy yy x y x x y y x y ∂ ∂ ∂ ∂ ∂ = − − − − − + ∂ ∂ ∂ ∂ ∂ ∂ U U U U U P MU A A A B B GU ɺɺ (2.21) where the vector P contains the consistent external tractions at the interfaces of the thin-layer (which result from the external tractions t and the body loads b ). The thin-layer matrices M , αβ A , α B and G are given by T 0 d h z ρ = ∫ M N N (2.22) T 0 d , , h z x y αα αα α = = ∫ A N D N (2.23) ( ) T 0 d h xy xy yx z = + ∫ A N D D N (2.24) T T 0 0 d d , , h h z z z z x y α α α α ′ ′ = − = ∫ ∫ B N D N N D N (2.25) T 0 d h zz z ′ ′ = ∫ G N D N (2.26) in which h is the thickness of the thin-layer and d d z ′ = N N . Appendix I tabulates the above matrices for an individual thin-layer consisting of a cross-anisotropic material and considering both a linear and quadratic interpolation, i.e., 2,3 m = , respectively. After the individual matrices are overlapped in the usual finite element sense (i.e. layer by layer and in the natural top down order of the interfaces), one obtains a narrowly banded set of global system matrices and vectors which characterizes the complete stack of thin-layers. The resulting system of partial differential equations has the same form as equation (2.21), but its shape is now block-tridiagonal and has a correspondingly larger number of equations. In the remaining part of the present chapter, equation (2.21) refers to the complete assembly of thin-layers. 2.3 Displacements in the wavenumber-frequency domain To solve the system of linear partial differential equations (2.21), the displacements U and tractions P are transformed from the space-time domain to the wavenumber-frequency domain by means of the triple Fourier transformations ( ) ( ) ( ) i , , , , e d d d x y t k x k y x y k k x y t x y t ω ω +∞ +∞ +∞ − − − −∞ −∞ −∞ = ∫∫∫ U U (2.27) ( ) ( ) ( ) i , , , , e d d d x y t k x k y x y k k x y t x y t ω ω ∞ ∞ ∞ − − − −∞ −∞ −∞ = ∫ ∫ ∫ P P (2.28) In the new domain, system (2.21) becomes ( ) ( ) 2 2 2 i x xx x y xy y yy x x y y k k k k k k ω   = + + + + + −   P A A A B B G M U (2.29) Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 23 where i 1 = − . All matrices in this expression are symmetric, except for x B and y B which are skew-symmetric. Although this system could be easily solved for U , it is both possible and convenient to first change the system of equations into a fully symmetric form by means of a similarity transformation. This is accomplished by multiplying every third row of the system (2.29) by i − and every third column by i . This operation solely affects the vectors P and U and the matrices x B and y B , leaving the other matrices unchanged. As a result of this transformation, the system of equations is now ( ) 2 2 2 x xx x y xy y yy x x y y k k k k k k ω   = + + + + + −   p A A A B B G M u ɶ ɶ ɶ ɶ (2.30) where p ɶ and u ɶ are obtained from P and U by multiplying every third row by i − . Also, x B ɶ and y B ɶ are obtained from x B and y B by reversing the sign of every third column [Note: in comparison with previous studies on the TLM (e.g. Kausel, 1986), in this work a reversed sign for the i factor is used for reasons of convenience]. After solving the system of equations (2.30), U is recovered by multiplying every third row of u ɶ by i and the displacements in the space-time domain can be obtained —at least formally— from the triple inverse Fourier transform ( ) ( ) ( ) ( ) i 3 1 , , , , e d d d 2 x y t k x k y x y x y x y t k k k k ω ω ω π ∞ ∞ ∞ − − −∞ −∞ −∞ = ∫ ∫ ∫ U U (2.31) The integrals in equation (2.31) can be evaluated numerically. However, by doing so the TLM loses its advantage over the integral transform techniques based on the transfer matrices or on the stiffness matrices, since these approaches, that also require the numerical evaluation of the inverse transformations, originate systems of equations of smaller size (usually, only the interfaces between the layers need to be discretized). In this way, a different procedure is followed, in which the displacement field is decomposed into a modal basis, similar to what is done in the modal superposition for linear dynamic analyses. As a result of this decomposition, the system of equations (2.30) can be diagonalized and that enables the evaluation in closed form expressions of at least one of the integrals of equation (2.31). If the integral to be evaluated is the outer integral, the fundamental solutions are obtained in the wavenumber-time domain (Kausel, 1994). If instead one changes the coordinates from Cartesian to cylindrical and evaluates the integral in the radial wavenumber, then the fundamental solutions are obtained in the space-frequency domain (Kausel and Peek, 1982; Kausel, 1981). Alternatively, if the inner integral of (2.31) is evaluated, then the fundamental solutions are obtained in a mixed space-wavenumber-frequency domain (2.5D domain), in which the plane-strain is the particular case 0 y k = (Barbosa and Kausel, 2012). In this work, in the propagation stage the geometry is assumed to be invariant and so the 2.5D fundamental solutions are of interest. On the other hand, in the reception stage the three-dimensionality of the problem has to be considered and so the cylindrical space-frequency domain solutions must be used. In the ensuing, the system of equations (2.30) is transformed in order to obtain the displacements U in the wavenumber-frequency domain through modal superposition. As a first step in that direction, the order of the degrees of freedom is rearranged, grouping first all horizontalx , then all horizontal- y and finally all verticalz degrees of freedom. This rearrangement is suggested solely to reveal the special structure possessed by the matrices in the system of equations (2.30) and the implications that the referred structures have on the Chapter 2 – Wave propagation in the soil: fundamental solutions 24 eigenvalue problems that are solved to find the modal basis. In practice, the degrees of freedom are ordered by interface and not by direction, which results in a reduction of the bandwidth of the matrices. Hence, after rearranging the degrees of freedom, the matrices ad vectors in system (2.30) attain the following structures x xx y z     =       A O O A O A O O O A y yy x z     =       A O O A O A O O O A x y xy x y −     = −       O A A O A A A O O O O O xz x T xz     =       O O B B O O O B O O ɶ y yz T yz     =       O O O B O O B O B O ɶ x y z     =       G O O G O G O O O G x y z     =       M O O M O M O O O M i x y z     =     −   u u u u ɶ i x y z     =     −   p p p p ɶ (2.32) where O is the null matrix, x y ≡ G G , x y z ≡ ≡ M M M , and xz yz ≡ B B . In addition, and except for the matrices xz B and yz B , all sub-matrices are symmetric and block-tridiagonal. Having rearranged the order of the degrees of freedom, it is now convenient to define the radial wavenumber k , the propagation angle ϑ , the transformation matrix T , its inverse 1 − T and the matrices A and C as 2 2 x y k k k = + cos x k k ϑ = sin y k k ϑ = (2.33) cos sin sin cos k ϑ ϑ ϑ ϑ     = −       I I O T I I O O O I 1 cos sin sin cos 1 k ϑ ϑ ϑ ϑ − −     =       I I O T I I O O O I (2.34) T x y xz z     =       A O O A O A O B O A 2 2 2 x x xz y y z z ω ω ω   −   = −     −   G M O B C O G M O O O G M (2.35) in which I is the identity matrix and O is a null matrix, both with the dimensions compatible with the submatrices in equation (2.32). By substituting each variable defined in equations (2.33)-(2.35) in the following equation, it can be shown that the system (2.30) is the same as ( ) -1 2 k = + p T A C Tu ɶ ɶ (2.36) or equivalently ( ) 2 k = + Tp A C Tu ɶ ɶ (2.37) The modal basis needed to decompose the displacements in a summation corresponds to the solution of the right eigenvalue problem in j k and j r Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 25 ( ) 2 j j k + = A C r 0 (2.38) Due to the structure of matrices A and C , the eigenvalue problem (2.38) can be decoupled into two eigenvalue problems, one in the x and z directions (generalized Rayleigh problem) and the other in the y direction (generalized Love problem): ( ) { } 2 2 T2 2 2 xj xx x xz Rj Rj zj xz z z z Lj y y y yj kk k ω ω ω       −     + =         −           + − = A O G M B 0 B A O G M A G M 0 f f f (2.39) In this way, the right eigenvalue problem (2.38) has two sets of eigenpairs: one set associated with the eigenvalues Rj k and right eigenvectors T T T Rj xj Rj zj k   =   r 0 f f and the other set associated with the eigenvalues Lj k and right eigenvectors T T Lj yj   =   r 0 0 f. Likewise, the left eigenvalue problem ( ) T 2 j j k + = l A C 0 (2.40) has two sets of eigenpairs: one set associated with the eigenvalues Rj k and left eigenvectors T T T Rj Rj xj zj k   =   l 0 f f and the other set associated with the eigenvalues Lj k and left eigenvectors T T Lj yj   =   l 0 0 f. The left and right eigenvectors satisfy the orthogonal conditions (Barbosa and Kausel, 2012) T Rj Rl jl Rj k δ =l Ar , T 3 Rj Rl jl Rj k δ = −l Cr T Lj Ll jl δ = l Ar , T 2 Lj Ll jl Lj k δ = −l Cr T 0 Rj Ll = l Ar , T 0 Rj Ll = l Cr (2.41) Having found the solutions of (2.38), the rotated displacements Tu ɶ are decomposed into a summation of the right eigenvectors, i.e. 1 1 R L N N Rj Rj Lj Lj j j = = = Γ + Γ ∑ ∑ Tu r r ɶ (2.42) where Rj Γ and Lj Γ are participation factors yet to be determined and R N and L N are the number of degrees of freedom in the Rayleigh and Love eigenvalue problems, respectively. After replacing the identity (2.42) in equation (2.37) and after pre-multiplying it by T Rl l , the latter becomes T T 2 1 1 R L N N Rl Rl Rj Rj Lj Lj j j k = =     = + Γ + Γ       ∑ ∑ l Tp l A C r r ɶ (2.43) Due to the orthogonal conditions expressed in (2.41), equation (2.43) is equivalent to ( ) T T 2 3 2 2 Rl Rl Rl Rl Rl Rl Rl Rl k k k k k k   = − Γ ⇔ Γ =   − l Tp l Tp ɶ ɶ (2.44) Chapter 2 – Wave propagation in the soil: fundamental solutions 26 If (2.37) is pre-multiplied instead by T Ll l , then T T 2 2 2 2 Ll Ll Ll Rl Ll Ll k k k k   = − Γ ⇔ Γ =   − l Tp l Tp ɶ ɶ (2.45) The combination of (2.44), (2.45) and (2.42) yields ( ) T T 2 2 2 2 1 1 R L N N Rj Lj Rj Lj j j Lj Rj Rj k k k k k = =       = +       − −     ∑ ∑ l Tp l Tp Tu r r ɶ ɶ ɶ (2.46) or equivalently ( ) T T 1 1 2 2 2 2 1 1 R L N N Rj Lj Rj Lj j j Lj Rj Rj k k k k k − − = =       = +       − −     ∑ ∑ l Tp l Tp u T r T r ɶ ɶ ɶ (2.47) Equation (2.47) can be further simplified into ( ) ( ) ( ) 2T T T 2 2 2 2 2 2 2 T T T 2 2 2 2 2 2 T 2 2 cos sin cos cos i sin cos sin sin i cos sin i i xj xj x xj xj y xj zj z xRj Rj Rj Rj xj xj x xj xj y xj zj z y Rj Rj Rj Rj Rj Rj zzj xj x Rj k k k k k k k k k k k k k k k k k k k k k ϑ ϑ ϑ ϑ ϑ ϑ ϑ ϑ ϑ + −   − − −       = + − − − −         + − p p p u p p pu up f f f f f f f f f f f f f f ( ) 1 T T 2 2 2 2 2T T 2 2 2 2 2 T T 2 2 2 2 1 1 sin sin cos sin cos cos R L N j zj xj y zj zj z Rj Rj yj yj x yj yj y Lj Lj N yj yj x yj yj y jLj Lj k k k k k k k k k k k k k ϑ ϑ ϑ ϑ ϑ ϑ ϑ = =           +         +− −       −   − −       − + − −           ∑ ∑ p p p p p p 0 f f f f f f f f f f f f (2.48) As a final step, taking into account the equality (Barbosa and Kausel, 2012) ( ) ( ) T T Rj zj xj zj xj Rj Rj Rj kk k k k k k k = − − f f f f (2.49) equation (2.48) can be written in the more convenient form Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 27 ( ) ( ) ( ) 2T T T 2 2 2 2 2 2 2 T T T 2 2 2 2 2 2 T 2 2 cos sin cos cos i sin cos sin sin i cos sin i i xj xj x xj xj y xj zj z xRj Rj Rj Rj xj xj x xj xj y xj zj zy Rj Rj Rj Rj zzj xj x Rj Rj k k k k k k k k k k k k k k k k k k k k k k ϑ ϑ ϑ ϑ ϑ ϑ ϑ ϑ ϑ ϑ + −   − − −       = + − − − −       +   − p p p u p p p u up f f f f f f f f f f f f f f ( ) 1 T T 2 2 2 2 2T T 2 2 2 2 2 T T 2 2 2 2 1 1 sin sin cos sin cos cos R L N j zj xj y zj zj z Rj Rj Rj yj yj x yj yj y Lj Lj N yj yj x yj yj y jLj Lj k k k k k k k k k k k k ϑ ϑ ϑ ϑ ϑ ϑ = =           +       +   − −       −   − −       − + − −           ∑ ∑ p p p p p p 0 f f f f f f f f f f f f (2.50) In the following, it will be implicitly understood that the eigenvalue problem for Rayleigh (shear vertical – pressure, SVP) waves will result in eigenvectors xj f , zj f whose components at the th m elevation and th j mode are written as ( ) ( ) , m m xj zj φ φ and their eigenvalues are j Rj k k =, while the eigenvectors yj f for Love (shear horizontal, SH) waves will have components written as ( ) m yj φ with eigenvalues j Lj k k =. In the light of equation (2.48), it is now convenient to define the set of kernels nj K given in Table 2.1. Table 2.1 : Kernels of fundamental solutions From equation (2.50) and in terms of the kernels in Table 2.1, the fundamental displacements ( ) ( ) , , mn x y U k k αβ ω at the th m elevation in direction α due to a unit load applied at the th n elevation in direction β can be expressed as listed in Table 2.2. ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 2 2 2 2 2 2 2 2 2 2 2 2 3 4 2 2 2 2 2 2 2 2 2 2 5 6 2 2 2 2 2 2 2 2 1 sin cos , cos sin , cos sin , x y j j j j j y x j j j j j j y x j j j j j j j j j j k k K K k k k k k k k k k K K k k k k k k k k k k k k k k K K k k k k k k k k k k k k ϑ ϑ ϑ ϑ ϑ ϑ = = = − − − = = = = − − − − = = = = − − − − Chapter 2 – Wave propagation in the soil: fundamental solutions 34 2.5.4 Vertical derivatives and internal stresses The z -derivatives can be obtained through the combination of the nodal displacements weighted by the derivatives of the associated shape functions. However, following that approach, the derivatives at the top and bottom interfaces of the thin-layers are not consistent with the tractions calculated in subsection 2.5.3, and so their accuracy is inferior. To compensate for the lack of precision of the vertical derivative, Kausel (2004) proposed an alternative strategy for their calculation. The procedure is based on the definition of secondary interpolation functions that are consistent with the tractions at the top and bottom interfaces of the thin-layer. In this work, that procedure is used to define the vertical derivatives and, subsequently, the internal stresses at the internal nodal interfaces. Consider the th i thin-layer (of expansion nn ), from which the displacements ( ) ( ) i j u αβ at the 1 nn + nodal interfaces, the horizontal derivatives ( ( ) ( ), i j x u αβ and ( ) ( ), i j y u αβ ), and the nodal tractions ( ( ) ( ) i j t αβ ) at the top and bottom interfaces are known for all directions , , x y z α = and for a source in the direction β . The tractions at the upper surface and the internal stresses at the same horizontal plane are related by T ( ) top top top (1) i xz yz zz β β β σ σ σ   =   t (2.73) while the tractions at the lower surface and the internal stresses at the corresponding plane are related by T ( ) bottom bottom bottom ( 1) i nn xz yz zz β β β σσσ +   = −   t (2.74) In their turn, the internal stresses and the derivatives of displacements are related by ( ) ( ) ( ) , , , , , , , xz t x z z x yz t y z z y zz t x x y y t z z G u u G u u u u D u β β β β β β β β β β σ σ σ λ = + = + = + + (2.75) The previous equation can be solved for the vertical derivatives, yielding ( ) ( ) ( ) , , , , , , , xz z x x z t yz z y y z t zz t x x y y z z t Gu uG Gu uG u u uD β β β β β β β β β β σ σ σ λ − = − = − + = (2.76) For each response direction α and for each source direction β , the values of the displacements at the 1 nn + nodal interfaces and of the vertical derivatives at the upper and lower interfaces are now known. In this way, it is possible to employ the 3 nn + known quantities and use Hermitian interpolation to define a polynomial of degree 2 nn + that approximates the vertical variation of the displacements. If the thin-layer is linear ( 1 nn = ) and its thickness is h , then Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 35 ( ) 2 3 u z A B z C z D z αβ αβ αβ αβ αβ = + + + (2.77) ( ) 2 , 2 3 z u z B C z D z αβ αβ αβ αβ = + + (2.78) 1( ) ( ) (2) (2) ( ) ( ) 2 3 (1) (1) ( ) ( ) 2 2 (2), (2), ( ) 2 3 3 2 2 (1), (1), 1 0 0 0 0 1 0 0 1 0 0 0 1 0 1 0 0 3 3 1 2 0 1 2 3 2 2 1 1 i i i i i i z z i z z Au u Bu u h h h Cu u h h h h Du u h h h h h h αβ αβ αβ αβ αβ αβ αβ αβ αβ αβ αβ αβ −                         = =         − − −         −             ( )i               (2.79) If instead the thin-layer is quadratic ( 2 nn = ), then ( ) 2 3 4 u z A B z C z D z E z αβ αβ αβ αβ αβ αβ = + + + + (2.80) ( ) 2 3 , 2 3 4 z u z B C z D z E z αβ αβ αβ αβ αβ = + + + (2.81) ( ) ( ) ( ) 1( ) (3) 2 3 4 ( ) (2) ( ) 2 3 4 (1) ( ) (3), ( ) 2 3 (1), 2 2 2 3 3 3 1 0 0 0 0 1 2 2 2 2 1 0 1 0 0 0 0 1 2 3 4 0 0 1 0 0 0 0 0 0 1 5 16 11 1 4 14 32 18 3 i i i i z i z Au Bu h h h h Cu h h h h Du Eu h h h h h h h h h h h h αβ αβ αβ αβ αβ αβ αβ αβ αβ αβ −                         = =                         − − − − − ( ) (3) ( ) (2) ( ) (1) ( ) 2 2 (3), ( ) 4 4 4 3 3 (1), 5 8 16 8 2 4 i i i i z i z u u u u h u h h h h h αβ αβ αβ αβ αβ                             − − −     (2.82) After finding the left-hand side of equations (2.79) and (2.82), the vertical derivatives of the displacements at the nodal interfaces can be determined using eq. (2.78) or eq. (2.81). As for the vertical derivatives ( ) ,z u z αβ of points in the interior of the thin-layer, though they can be calculated using the same two equations, one chooses to use instead the original interpolation functions to determine these variables, i.e., ( ) ( ) 1( ) , ( ), 1 nn i z j j z j u z N z u αβ αβ + = = ∑ (2.83) With all the first derivatives known, eq. (2.6) can be used to calculate the internal stresses. Regarding the calculation of the second derivatives involving the z direction ( 2 x z ∂ ∂ ∂ , 2 y z ∂ ∂ ∂ and 2 2 z ∂ ∂ ), , yz u αβ is obtained by multiplying , z u αβ by i y k − , while , xz u αβ is obtained by differentiating equation (2.76) with respect to x , i.e. Chapter 2 – Wave propagation in the soil: fundamental solutions 36 ( ) ( ) ( ) , , , , , , , , , , xz x z xx x xz t yz x z xy y xz t zz x t x xx y xy z xz t Gu uG Gu uG u u uD β β β β β β β β β β σ σ σ λ − = − = − + = (2.84) The traction derivatives , z x α β σ are calculated as explained in subsection 2.5.3, with the exception that matrices ( ) p jR β Λ and ( ) p jL β Λ must be replaced with matrices ( ) 1 i p jR β + − Λ and ( ) 1 i p jL β + − Λ , respectively. For the calculation of the matrices (3) jR β Λ and (3) jL β Λ the integrals (3) nj I of the form ( ) i (3) 3 1 2 e d , 1,...,6 x k x nj x nj x I x k K k n π +∞ − −∞ = = ∫ (2.85) are needed. Closed form expressions for these integrals are given in Table 2.6. Table 2.6: Closed form expressions for (3) nj I ( 2 2 Im 0 j y k k − < ) ( ) ( ) { } ( ) ( ) { } 2 2 2 2 2 2 2 2 i i (3) 3 1 1 3 2 3 i i (3) 3 2 2 2 2 2 2i i (3) 3 2 2 4 3 3 2 (3) 4 1 2 1 2 1 2 1 2 e d sgn e 2i e d e i e 2i sgn e d e e 2i j y x j y y x j y y x x k k j y k x j x j x x k k k x y k x j x j x j y y j x k k k x k x j x j x j y y j j k k I k K k x k I k K k k k k k x I k K k k k k k I π π π +∞ − − − −∞ +∞ − − − − −∞ +∞ − − − − −∞ − = = = = − − = = − − = ∫ ∫ ∫ ( ) ( ) { } ( ) ( ) ( ) 2 2 2 2 2 2 i i 3 2 2 2 4 42 3 2 2 2 i i (3) 3 5 5 2 2 i i (3) 3 6 6 1 2 1 2 sgn e d e e 2i e d e 2i e d sgn e 2i j y y x j y x j y x x k k k x k x x j x y j y y j j y x k k k x j x j x j y j y x k k k x j x j x j x k K k k k k k k k k I k K k k k k k I k K k x k π π π +∞ − − − − −∞ +∞ − − − −∞ +∞ − − − −∞ = − + − = = − = = ∫ ∫ ∫ [Note about Table 2.5 and Table 2.6: ( ) 3 2 2 2 j y k k− must be calculated as ( ) 2 2 2 2 j y j y k k k k − − , with 2 2 Im 0 j y k k − < ; a direct use of the expression ( ) 3 2 2 2 j y k k− might possibly assign the wrong sign to the result.] The remaining second derivative, , zz u αβ , can be calculated resorting to the Navier equation (Achenbach, 1973) ( ) 2 , , i jj j ij i i G u G u F u λ ρω + + + = − (2.86) which after being solved for , zz u αβ yields (it is assumed that 0 i F = , i.e., no internal sources) Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 37 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 , , , , , , 2 , , , , , , 2 , , , , , 2 x x xx y xy z xz x xx x yy x zz y x xy y yy z yz y xx y yy y zz z x xz y yz z xx z yy z zz u G u u u G u u uG u G u u u G u u uG u G u u G u u u G β β β β β β β β β β β β β β β β β β β β ρω λ ρω λ ρω λ λ − − + + + − + = − − + + + − + = − − + + − + = + (2.87) 2.5.5 Validation To validate the equations derived in this and previous sections, the response of a homogeneous full-space obtained with the TLM is compared with the analytical solution derived by Tadeu and Kausel (2000). The isotropic full-space has mass density 1 ρ = , shear modulus 1 G = and Poisson’s ratio 0.25 ν = , and is subjected to a time harmonic line load of the form ( ) ( ) ( ) ( ) , , , exp i i y x y z t x z t k y δ δ ω = −b, being the excitation frequency 1 f = Hz ( 2 ω π = ). The thin-layer model used to simulate the full-space consists of an elastic layer with thickness 20 s λ divided into 200 thin-layers of quadratic expansion — a discussion on discretization errors can be found in Park and Kausel (2004b) — and supplemented at its upper and lower horizons with paraxial boundaries (Seale and Kausel, 1989), which are used to simulate the infinite domain — in section 2.6 of this chapter, a more efficient strategy to model infinite domains is presented. The parameters s λ is the wavelength of the shear wave, which is defined by s s s C G C f λ ρ = = (2.88) The load is applied at the middle surface of the elastic layer. To avoid strong oscillations in the response, a small amount of damping is considered ( 0.005 P S ξ ξ = = ), which renders the wave velocities complex. In the first validation scenario, the displacements induced by a vertical load are computed as function of the wavenumber y k at the horizontal plane 0 z = and horizontal distances 0.01 s x λ =, s x λ = and 5 s x λ =. Figure 2.2 depicts the comparison between the vertical displacements obtained with the TLM and those calculated with the analytical solution. As can be observed in Figure 2.2b) and c), the match between the exact solution and the TLM solution is very good. Nonetheless, very close to the load (Figure 2.2a), one can observe a rather small difference in the real part, which is due to discretization effects. This is because the thickness of the thin-layers is only 0.1 s λ while the receiver is placed at one tenth of that distance from the source. Still, given the excellent quality of the comparison even at that short range, the results obtained clearly demonstrate the robustness of the TLM solution. Observe that at large distances the response decays very fast with the wavenumber y k beyond the threshold / s s k C ω = (the branch point), while below that value the response is highly wavy. Hence, when computing the inverse transform from y k –space into y –space for remote points, one can truncate the integrals at the branch point, but then again because of the rapid oscillations one must consider a sufficiently dense number of points below that threshold. Conversely, for receivers at close range, the response functions are less wavy, but they also Chapter 2 – Wave propagation in the soil: fundamental solutions 38 decay more slowly with y k . Hence, their Fourier inversion must include points beyond the branch point even if one can get away with coarser spacing. Figure 2.2: Vertical displacements at a) 0.01 s x λ =, b) s x λ = , and c) 5 s x λ =. Solid lines = TLM solution (real part — blue; imaginary part — red). Circles = analytical solution In the second validation example, both displacements and derivatives are computed as function of the horizontal distance x and considering all three directions for the load and for the response. The longitudinal wavenumber is 0.4 y k= and the receivers are placed at the depth 4 s z λ = and up to the horizontal distance max 4 s x λ =. The stresses are not compared because they can be calculated based on the derivatives of displacements, and if the latter are correct, then the former are also correct. Likewise, the y -derivatives are not represented because they result from the multiplication of other response fields (displacement or derivative) by i y k − . Figures 2.3 to 2.8 depict the comparison between the theoretical solution and the responses obtained with the TLM. All figures suggest that the two approaches yield results that are virtually the same, thus confirming that the TLM is indeed capable of reproducing the wave motion with high accuracy. This validates the expressions derived in the previous sections of this work. 0 2 4 6 8 10 -0.5 0 0.5 1 k y u zz a) 0246810 -0.2 -0.1 0 0.1 0.2 k y u zz b) 0 2 4 6 8 10 -0.04 -0.02 0 0.02 0.04 k y u zz c) Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 39 Figure 2.3: Displacements u αβ . Solid lines = TLM solution (real part — blue; imaginary part — red). Circles = analytical solution Figure 2.4: Displacement derivatives , x u αβ . Solid lines = TLM solution (real part — blue; imaginary part — red). Circles = analytical solution 024 -0.05 0 0.05 u xx 0 2 4 -4 -2 0 2x 10 -3 u xy 0 2 4 -0.02 0 0.02 u xz 024 -4 -2 0 2x 10 -3 u yx 02 4 -0.05 0 0.05 u yy 0 2 4 -5 0 5x 10 - 3 u yz 024 -0.02 0 0.02 u zx 024 -5 0 5x 10 -3 u zy 02 4 -0.02 0 0.02 u zz 024 -0.1 0 0.1 u xx,x 0 2 4 -10 -5 0 5x 10 -3 u xy,x 0 2 4 -0.1 0 0.1 u xz,x 024 -10 -5 0 5x 10 -3 u yx,x 0 2 4 - 0.2 0 0.2 u yy,x 0 2 4 -5 0 5 10 x 10 -3 u yz,x 024 -0.1 0 0.1 u zx,x 0 2 4 -5 0 5 10 x 10 -3 u zy,x 0 2 4 -0.1 0 0.1 u zz,x Chapter 2 – Wave propagation in the soil: fundamental solutions 40 Figure 2.5: Displacement derivatives , z u αβ . Solid lines = TLM solution (real part — blue; imaginary part — red). Circles = analytical solution Figure 2.6: Displacement derivatives , xx u αβ . Solid lines = TLM solution (real part — blue; imaginary part — red). Circles = analytical solution 024 -0.5 0 0.5 u xx,z 0 2 4 -5 0 5 10 x 10 -3 u xy,z 0 2 4 -0.1 0 0.1 u xz,z 024 -5 0 5 10 x 10 -3 u yx,z 0 2 4 -0.5 0 0.5 u yy,z 0 2 4 -0.01 0 0.01 0.02 u yz,z 024 -0.1 0 0.1 u zx,z 0 2 4 -0.01 0 0.01 0.02 u zy,z 0 2 4 -0.1 0 0.1 u zz,z 024 -0.5 0 0.5 u xx,xx 0 2 4 -0.05 0 0.05 u xy,xx 0 2 4 -0.5 0 0.5 u xz,xx 024 -0.05 0 0.05 u yx,xx 0 2 4 -1 0 1 u yy,xx 0 2 4 -0.05 0 0.05 u yz,xx 024 -0.5 0 0.5 u zx,xx 0 2 4 -0.05 0 0.05 u zy,xx 0 2 4 -0.5 0 0.5 u zz,xx Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 41 Figure 2.7: Displacement derivatives , xz u αβ . Solid lines = TLM solution (real part — blue; imaginary part — red). Circles = analytical solution Figure 2.8: Displacement derivatives , zz u αβ . Solid lines = TLM solution (real part — blue; imaginary part — red). Circles = analytical solution 024 -0.5 0 0.5 u xx,xz 0 2 4 -0.05 0 0.05 u xy,xz 0 2 4 -0.5 0 0.5 u xz,xz 024 -0.05 0 0.05 u yx,xz 0 2 4 -1 0 1 u yy,xz 0 2 4 -0.05 0 0.05 u yz,xz 024 -0.5 0 0.5 u zx,xz 0 2 4 -0.05 0 0.05 u zy,xz 0 2 4 -0.5 0 0.5 u zz,xz 024 -1 0 1 2 u xx,zz 0 2 4 -0.05 0 0.05 u xy,zz 0 2 4 -0.5 0 0.5 u xz,zz 024 -0.05 0 0.05 u yx,zz 0 2 4 -1 0 1 2 u yy,zz 0 2 4 -0.1 0 0.1 u yz,zz 024 -0.5 0 0.5 u zx,zz 0 2 4 - 0.1 0 0.1 u zy,zz 0 2 4 -0.5 0 0.5 u zz,zz Chapter 2 – Wave propagation in the soil: fundamental solutions 42 2.6 Modeling unbounded domains The TLM is a semi-discrete numerical technique for the analysis of wave motion in layered media and relies on a finite element discretization in the direction of layering. Due to the discrete character of the TLM, by itself the analyses are limited to bounded domains. Nevertheless, in the mid eighties, the Paraxial Boundaries (PB) were coupled to the TLM to allow the simulation of infinite domains (Seale and Kausel, 1989). Very briefly, the PB for the TLM can be obtained by expanding the stiffness matrices (Kausel and Roesset, 1981) in a Taylor series in the wavenumber and retaining only the first three terms. Hence, this technique works very well for small wavenumbers and not so well for higher wavenumbers. In other words, waves propagating vertically or almost vertically are mostly absorbed when they reach the paraxial boundary while waves propagating with a considerable horizontal component are mostly reflected. For this reason, the PBs are usually augmented with buffer (elastic) layers that are thick enough so that the component of waves that is reflected at the PB returns to the region of interest at a horizontal coordinate larger than the maximum distance of interest. Explicit expressions for the PB matrices can be found in the thesis of Park (Park, 2002, p. 284 for SH waves and p. 289 for the SVP waves) while comments and considerations about their stability can be found in the works by Kausel (1988, 1992). More recently, Barbosa et al. (2012) successfully coupled the Perfectly Matched Layer (PML) to the TLM, resulting in a more efficient technique to model unbounded domains. The PML is a numerical technique used for purposes similar to those of absorbing or transmitting boundaries, namely to suppress undesirable echoes and reflections of waves in infinite media modeled with discrete finite systems. The technique was introduced in the nineties by Berenger (1994), who developed and coupled it to the time-domain finite differences method for the analysis of electromagnetic fields. Initially, the formulation of the PML followed the split-field approach, in which the systems of equations are solved both for displacements and stresses, but later new formulations for the PML were derived, namely by considering the PML as an equivalent anisotropic material (Gedney, 1996; Teixeira and Chew, 1998) or by stretching the coordinates to the complex space (Chew and Weedon, 1994; Hugonin and Lalanne, 2005). The PML has also been applied to elastodynamic problems, both in time and frequency domains (Basu and Chopra, 2003; Basu and Chopra, 2004; Basu, 2009; Harari and Albocher, 2006). A good literature review on the subject can be found in (Kucukcoban and Kallivokas, 2010). In the work (Barbosa et al., 2012), the PML is formulated based on the coordinate stretching approach. This approach consists in stretching the real space to a complex space by means of position-dependent complex-valued scaling functions, which begin with unit values at the interface or horizon delimiting the elastic region and then attain progressively larger complex values with the distance from this horizon, which causes the waves within the PML to attenuate exponentially (Johnson, 2008). Since there is no impedance contrast at the PML boundary, no reflections take place no matter what the angle of propagation of the waves entering the PML is. In the following subsections, it is shown how the coordinate stretching allows the simulation of infinite domains and then the PML is coupled to the TLM. 2.6.1 Coordinate stretching Consider a plane wave travelling at an angle θ with the vertical direction ( z ) in a medium whose wave speed is C . This wave has the form Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 43 ( ) i sin cos , , e t x z C C u x z t A ω ω ω θ θ   − −     = (2.89) No restriction is made regarding the vertical coordinate, and thus admit that this coordinate may assume the complex values, and denote it by z . Equation (2.89) is still valid but now z must be replaced by z . Assume also that the imaginary part of z depends on the depth (real part) and define z as ( ) i z z z = − Ψ (2.90) where ( ) z Ψ is a function yet to be determined. After replacing equation (2.90) into (2.89), the latter becomes ( ) ( ) i sin cos cos , , e e t x z z C C C u x z t A ω ω ω ω θ θ θ   − − −Ψ     = (2.91) The aim is to attenuate the waves that enter a finite PML region defined by 0 H z > > (Figure 2.9). For waves that propagate in the positive z direction ( pos cos 0 θ > ), for the amplitude of the waves to decay as z increases, the exponential term ( ) ( ) exp cos z C ω θ − Ψ must decrease and consequently ( ) z Ψ must increase with z . Similarly, for waves that propagate in the negative z direction ( neg cos 0 θ < ), for their amplitude to decay in the direction of propagation, ( ) z Ψ must obey the same rule. A possible choice for ( ) z Ψ that respects the established rule is ( ) ( ) 0 z z d ψ ζ ζ Ψ = ∫ (2.92) where ( ) z ψ is an always positive stretching function. Figure 2.9: Propagation of a wave inside the PML region In theory, ( ) z ψ might assume any shape as long as ( ) ( ) 0 H Ψ < Ψ (Bienstman and Baets, 2002). However, once the domain is discretized so that the differential equations can be solved, spurious reflections occur due to changes in ( ) z ψ and therefore it is convenient that this function changes smoothly with z . A commonly used stretching function is (Basu and Chopra, 2003) ( ) 0 m z zH ω ψω   =     (2.93) H z neg pos θ π θ = − pos θ Chapter 2 – Wave propagation in the soil: fundamental solutions 50 Bienstman et al., 2001), but in these works, only one branch exists. In the present work, the two branches are justified by the existence of two different body waves. The branches start at the wavenumbers s s k C ω = and p p k C ω = and the number of modes contained in each branch equals the number of degrees of freedom associated with the PML layer. Vertical displacements due to a vertical line load ( 0 y k = ) The modes obtained in subsection 2.6.3.1 are now combined in order to obtain the displacements in the interface between the layer and the half-space due to loads at the same elevation. Figures 2.14 and 2.15 plot the displacements for the frequencies 0.2 Hz f = and 1Hz f = obtained both with the TLM and with numerical integration on the wavenumbers of the displacements given by the stiffness matrices. Figure 2.14: Displacements of the layered half-space: TLM vs Stiffness matrices ( 0.2 Hz f = ) Figure 2.15: Displacements of the layered half-space: TLM vs Stiffness matrices ( 1Hz f = ) 0 5 10 15 20 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Distance to the source u zz Re(u zz ) - TLM Im(u zz ) - TLM Re(u zz ) - Stiff. Mat. Im(u zz ) - Stiff. Mat. 0 5 10 15 20 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 Distance to the source u zz Re(u zz ) - TLM Im(u zz ) - TLM Re(u zz ) - Stiff. Mat. Im(u zz ) - Stiff. Mat. Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 51 The comparison of the results obtained with the two approaches shows that except for very small shifts, the displacements agree very well. Hence, it is concluded that the PML is accurate in reproducing the behavior of infinitely deep stratified domains. 2.7 Solution of the eigenvalue problems One of the more time consuming and complex steps in the TLM process is the solution of the eigenvalue problems (2.39). Even though commercial software, like Matlab, contain functions capable of solving both the linear and the quadratic eigenvalue problems with complex matrices (for example, eig and polyeig), these functions do not take advantage of the banded structure of the matrices and therefore become inefficient for the solution of eigenproblems with dimension of just a few hundred. Also, these routines, which are based on the QZ algorithm (Kressner, 2005) and thus based on iterative rotations of the modal basis until convergence, may yield eigenvectors that do not respect entirely the orthogonal conditions (2.41) due to accumulation of round-of errors, hence introducing errors in the remaining steps of the modal combinations: though rare, this case has been observed when a large number of thin-layers and PMLs are used. Alternative approaches to the QZ algorithm are the methods based on the Power Method. This family of methods is based on successive matrix-vector multiplication until the direction of the resulting vector has converged to the eigen direction. Several variations exist, namely the Power Method itself, the Inverse Iteration method, in which instead of the matrix-vector multiplication, a system of equations is solved, and the Inverse Iteration with Rayleigh Shift method (Shit and Invert method), in which the eigenvalues are shifted in order to speed up convergence (Saad, 1992). In this family of methods, the eigenvectors are found one by one, being then deflated from the eigen base (remove them from the modal spectrum) in order to avoid convergence to the same pair. Projection methods can also be employed: the most common ones are the Lanczos method (or Arnoldi method, for non-hermitian matrices) and the Locally Optimal Preconditionated Conjugated Gradients method. These projection methods iterate with a group of vectors instead of a single vector and are used to extract a small spectrum range, not the full spectrum of eigenpairs (Saad, 1992). Because the complete spectrum of the system of matrices is needed, the Inverse Iteration with Rayleigh Shift is chosen for the solution of the eigenvalue problems (2.39). In the next subsections, it is explained how to apply this method taking into account the special structure of the TLM matrices. 2.7.1 Eigenvalue problem for SH waves The linear eigenvalue problem to be solved has the form ( ) 2 j j k + = A C v 0 (2.109) where matrices A and C are complex, symmetric and banded. The eigenvalues j k and eigenvector j v are also complex. The bandwidth of the matrices is 2 if the thin-layers are of linear expansion and is 3 if the thin-layers are of quadratic expansion. The eigenvalue problem admits 2 N solutions, being N the dimension of the matrices, and if the pair ( ) , j j k v is a solution of (2.109), then the pair ( ) , j j k− v is also a solution of (2.109). In this way, only Chapter 2 – Wave propagation in the soil: fundamental solutions 52 half of the spectrum of the system needs to be computed. Taking into account the orthogonal conditions expressed in the second row of eq. (2.41), the Inverse Iteration with Rayleigh Shift method applied to the linear pencil (2.109) results in the algorithm described in Table 2.7. Table 2.7: Algorithm for the solution of the linear eigenvalue problem ( ) rand j N = v Initial guess T j j j i = − v v v Av Normalize the initial guess against previously found eigenvectors ( 1,2,... 1 i j = − ) j j = w Av Auxiliary vector T T j j j j j k= − v Cv v w Rayleigh quotient Iterate (until Rayleigh quotient has converged) 2 j k = + E A C From shifted matrix 1 j j − = − v E w New approximation of eigenvector T j j j i = − v v v Av Normalize against previously found eigenvectors (if j i k k ≈ , 1,2,... 1 i j = − ) T 2 T j j j j j j k k= + v w v Av Update the Rayleigh quotient j j = w Av Auxiliary vector T j j j j j j δ δ δ =  =  =   v w v v w w Normalize the approximation The normalization step (against previously found eigenpairs) is forced at the beginning of the procedure and during the iterative procedure in order to avoid finding the same solution twice. Notice however that at the iteration loop, the normalization is only performed against the eigenvectors whose associated eigenvalue approximates the current Rayleigh quotient. Another important aspect of the algorithm is the solution of the system of equations 1 j j − = − v E w . Because matrices A and C are symmetric and banded, matrix E contains similar structure. Hence, the application of Gaussian elimination to this system of equations is very efficient. In opposition to the Power method or the Inverse Iteration method, the chosen approach requires the calculation and factorization of matrix E at every iteration. Nevertheless, because the cost of factorization of the matrix is very low (due to its slim banded structure), because the convergence is greatly improved (due to the application of the Rayleigh shift), and Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 53 because the normalization (deflation) is imposed only against the pairs that are in close proximity, then the Inverse Iteration with Rayleigh Shift method turns out to be the most efficient tool among the three mentioned methods. One last comment: the speed of convergence depends on the quality of the initial guess of j v . In this work, the initial guess is a random value. However, because the eigenvalue problem needs to be solved for different values of ω ( 2 ω = − C G M ), if two successive frequencies are close enough to each other ( 2 1 d ω ω ω = + ), one can use the eigenvectors computed for 1 ω as initial guesses for the eigenvectors of 2 ω . This strategy shall improve the convergence of the procedure (examples have shown that the average number of iteration per eigenpair reduces from 7-8 to 2-3; this strategy has not been used in this work due to complications that occur when two eigenvalues are too close from each other, and as a consequence the solutions jump from one branch to another). 2.7.2 Eigenvalue problem for SVP waves The first eigenvalue problem of eq. (2.39) is equivalent to ( ) 2 j j j k k + + = A B C v 0 (2.110) where x z   =     A Ο A Ο A xz T xz   =     ΟB BB Ο 2 2 x x z z ω ω   − =   −   G M Ο CΟG M xj v zj v   =     f f Rj j k k = Taking into consideration the original order of the degrees of freedom instead of the order organized by direction (see section 2.3), matrices A , B and C result symmetric and banded, being the bandwidth 4 if the thin-layers are of linear expansion or 6 if the thin-layers are of quadratic expansion. System (2.110) can be further expanded and written as j j j j j j j k k k−         =                 v v B C A Ο v v CΟ Ο C (2.111) and so it can be easily concluded that the eigenvalue problem has 2 N solutions, being N the dimension of matrices A , B and C . Also, due to the symmetric properties of these matrices and due to the special structure of matrix B , if the pair ( ) , j j k v is a solution of (2.110), then the pair ( ) , j j k− v is also a solution. The modified eigenvector j v coincides with j v , except for the components associated with the z direction, which have their sign reversed. Because the quadratic eigenvalue problem is slightly more complex than the linear one, some explanations are needed before presenting the solution procedure. Of equation (2.111), drop the modal index j and consider a general iteration (say, the th l iteration) of the Inverse Iteration with Rayleigh Shift: Chapter 2 – Wave propagation in the soil: fundamental solutions 54 1 1 l l l l l l k k + + +−         =         −        A B C u u AΟ C C v v ΟC (2.112) Handling simultaneously both rows of the system (2.112), it is possible to obtain ( ) 1 1 2 1 l l l l l l l l l l k k k k + + + − =   + + = −  u v v A B C u Cv Au (2.113) Notice that at each iteration, it is needed to solve a system of equations with dimension N and not 2 N , as could be suggested by equation (2.112). As for the Rayleigh quotient, its value is given by T T T T T T T T 2 l l l l l l l l l l l l l l l l l k             +     = = −−                 u B C u v v C O u Bu u Cv u A O v Cv u Au u v v O C (2.114) However, in this work, instead of calculating the Rayleigh quotient according to (2.114), the quotient is updated in every iteration according to the recursive formula (Waas, 1972) T T 1 1 1T T 1 1 1 1 l l l l l l l l l l k k + + + + + + + − = + − u Au v Cv u Au v Cv (2.115) Equation (2.115) works as well as eq. (2.114), but is more convenient in terms of computational resources. The orthogonality between a iteration of the th j eigenvector and an eigenvector already converged (the th i eigenvector) is imposed through ( ) ( ) T T T T j j i i i i j i i j i j i i j i j j i i k k k k −         = − − − +                 u u v v v Cv u Av v Cv u Av v v v v (2.116) As can be seen in eq. (2.116), both the eigenvector associated with i k and the eigenvector associated with i k − are being deflated. The algorithm used for the solution of the eigenvalue problem takes into account the major steps described above and is listed in Table 2.8. Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 55 Table 2.8: Algorithm for the solution of the quadratic eigenvalue problem ( ) rand j N =v Initial guess ( ) rand j N =u Initial guess (auxiliary) Equation (2.116) Normalize the initial guess against previously found eigenvectors ( 1,2,... 1 i j = − ) j j = w Au Auxiliary vector j j = y Cv Auxiliary vector T T T T 2 j j j j j j j j l k+ =− u Bu u y v y u w Rayleigh quotient Iterate (until Rayleigh quotient has converged) 2 j j k k = + + E A B C From shifted matrix ( ) 1 j j j j k − = − u E y w ( ) j j j j k = −v u v New approximation of eigenvector Equation (2.116) Normalize against previously found eigenvectors (if j i k k ≈ or j i k k ≈ − , 1,2,..., 1 i j = − ) T T T T j j j j j j j j j j k k − = + − u w v y u Au v Cv Update the Rayleigh quotient j j = w Au Auxiliary vector j j = y Cv Auxiliary vector ( ) T 2 2 T j j j j j j j j j j j k δ δ δ δ − =  =  =  =  u w u Cu v v w w y y Normalize the approximations 2.8 Conclusions In this chapter, the Thin-Layer Method is introduced as a tool to calculate the response of layered domains (for example, soil) to dynamic loads. The TLM is extended to the 2.5D domain and closed form expressions for the fundamental displacements, derivatives, tractions and stresses are given and validated. It is worth noting that the proposed methodology relies on the solution of two eigenvalue problems that in no way depend on the horizontal wavenumbers x k and y k , hence they need to be solved only once for each frequency. The 3D Chapter 2 – Wave propagation in the soil: fundamental solutions 56 fundamental solutions obtained with the TLM can be found in other works on the subject (Kausel, 1981). A new and more efficient procedure, the Perfectly Matched Layer, is proposed for the simulation of infinite domains. The PML is used to solve examples of full-spaces and layered half-spaces and is validated by means of the same examples. This new procedure is much more efficient than the previous one (the paraxial boundaries), and from now on it is recommended that the PB be replaced by the PML. In the last section of this chapter, efficient algorithms for the solution of the two eigenvalue problems are described. Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 57 3. Numerical tools for soilstructure interaction 3.1 Introduction Soil-structure interaction received a lot of attention during the XX th century. It is an interesting and complex problem, and many works covering a wide range of fields can be found in the literature. The interested reader is referred to the work (Kausel, 2010) for a historical overview on the subject. For the case of traffic induced vibrations, track and building are coupled through the ground, and so, both in the propagation stage and in the reception stage, the interaction between the soil and the track or the building must be considered. Due to the differences in the typology, the ideal numerical tool to model each of the sub-domains may vary. For example, to model the soil, which is an infinite domain, the Boundary Element Method (BEM) is preferred because it can account for the radiation of waves to the infinity (Dominguez, 1993). On the other hand, the use of the Finite Element Method (FEM) reveals to be more appropriate to model the behavior of buildings and tracks, since these structures are circumscribed and generally irregular. Additionally, in the propagation stage the geometry of the problem can be assumed invariant in the longitudinal direction, and thus a 2.5D formulation is advantageous, while in the reception stage the problem is limited to the structure under analysis and the surrounding environment, and therefore a 3D formulation is more appropriate. In this work, the propagation stage and the reception stage are decoupled. Hence, the wave field that the track induces in the soil is calculated disregarding the presence of the building in the far field. This simplification has been used by other authors (François, 2008; Fiala et al., 2007) and is valid when the distance between the building and the track is larger than the characteristic wavelengths of the soil, i.e., the simplification is valid for the medium and the high frequency range. Nevertheless, in this work this assumption is also used for the low frequencies. In the present chapter, the 3D BEM, the 2.5D BEM and the 2.5 D FEM are described. 3.2 3D Boundary Element Method 3.2.1 Introduction The Boundary Element Method is a discrete numerical procedure that can be used to solve partial differential equations. This procedure relies on the discretization of the boundary of the domain, as opposed to the Finite Element Method, for which the whole domain must be modeled. For this reason, for problems involving very large domains (such as the soil) the BEM becomes advantageous over the FEM as it avoids the discretization of the interior of the domain, which results in a reduced number of degrees of freedom. Another advantage of the BEM is that it deals with the radiation of waves to infinity exactly, contrarily to the FEM, in which special procedures need to be applied at the boundaries of the truncated domain in order to model unbounded domains. As drawbacks, the BEM renders systems of equations which involve fully populated and non-symmetric matrices and requires the knowledge of the Chapter 3 – Numerical tools for soil-structure interaction 58 Fundamental Solutions of auxiliary domains, which can be homogeneous full-spaces, halfspaces or layered spaces. In this work the soil is assumed to be horizontally stratified — in the work of Jones (2010), the importance of the inclination of layers is assessed. For this reason, the more appropriate auxiliary domain is the layered half-space. Using the fundamental solutions of such auxiliary domain, the BEM requires solely the discretization of the surfaces of the soil interacting with the structures. Unfortunately, such fundamental solutions are not known in closed form expressions. An alternative would be to use the full-space fundamental solutions, for which analytical expressions are known. However, when considering such auxiliary domain, the free-surface of the soil and of the interfaces between the different layers that characterize the soil must also be discretized, and that results in a substantial increase of the number of degrees of freedom, which renders the approach less attractive. Three-dimensional formulations of the BEM can be found in the works of Brebbia and Dominguez (1992) for elastostatic problems and of Dominguez (1993) for elastodynamic problems. For train induced vibration problems, the time domain 3D BEM has been used by Galvín (2007) and O'Brien and Rizos (2005), who considered the fundamental solutions of homogeneous full-spaces, by François (2008), who considered the fundamental solutions of layered half-spaces, and by Bode et al. (2002), who considered the fundamental solutions of homogeneous half-spaces. The frequency domain 3D BEM has been used by Hubert et al. (2001), who considered full-space fundamental solutions, and by Fiala et al. (2007), who considered the layered half-space fundamental solutions. In this work, a frequency domain 3D BEM procedure is coupled to the 3D FEM in order to obtain the response of a structure due to an incoming wave field. In the following section, the 3D BEM is presented. 3.2.2 Integral representation Consider an elastic body Ω with boundary Γ and two elastodynamic states described by the displacement fields ( ) , k u t x and ( ) * , k u t x , the initial displacements ( ) 0k u x and ( ) * 0k u x , the initial velocities ( ) 0k v x and ( ) * 0k v x , the body forces ( ) , k b t ρ x and ( ) * , k b t ρ x and the boundary tractions ( ) , k p t x and ( ) * , k p t x , where 1, 2, 3 k = corresponds to the x , y and z directions and where x is a point with coordinates ( ) , , x y z =x . The two elastodynamic states are interrelated by the Reciprocal Theorem , which states (Achenbach, 1973) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) { } ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) { } * * * * 0 0 * * * * 0 0 , , d , , , , d , , d , , , , d k k k k k k k k k k k k k k k k p t u t b t u t u u t v u t p t u t b t u t u u t v u t ρ ρ Γ Ω Γ Ω ∗ Γ + ∗ + + Ω = = ∗ Γ + ∗ + + Ω ∫ ∫ ∫ ∫ x x x x x x x x x x x x x x x x ɺ ɺ (3.1) In equation (3.1), a dot over a variable represents the derivative with respect to time and the operator “ * ” represents the time convolution, defined by ( ) ( ) ( ) ( ) ( ) ( ) 0 0 * d d t t f t g t f t g g t f τ τ τ τ τ τ = − = − ∫ ∫ (3.2) Also, the Einstein notation is used, i.e., a repeated index in a term implies the summation over that index. For example, 3 * * * * * 1 1 2 2 3 3 1 k k k k k p u p u p u p u p u = ∗ = ∗ = ∗ + ∗ + ∗ ∑ (3.3) Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 59 Now, assume that the second elastodynamic state (the one associated with the supperscript * ) is elicited by a body force of the form ( ) ( ) * , ,0 k kl b t ρ δ δ = x ξ (3.4) where 1, 2, 3 l = is the direction of the body force, kl δ is the Kronecker delta, ( ) δ … is the Dirac delta function and ( ) , , x y z ξ ξ ξ =ξ is the point where the body force is applied, which can be inside or outside the domain Ω . In this case, the displacements induced by such load correspond to the fundamental displacements in the time domain and, in the ensuing, are denoted by ( ) * , , kl u t x ξ , in which the first index represents the direction of the displacement and the second index represents the direction of the body load as given in equation (3.4). Analogously, the fundamental tractions at the boundary Γ are denoted by ( ) * , , kl p t x ξ and are calculated with ( ) ( ) 3 * * 1 , , , , kl kjl j j p t t n σ = = ∑ x ξ x ξ (3.5) where j n are the components of the vector n that is normal to the boundary Γ at the point x (pointing outwards) and where ( ) * , , kjl t σ x ξ are the fundamental stresses in the time domain induced at the point x by a point load at ξ . Assume also that the first elastodynamic state presents no body forces and that it is initially at rest, i.e., ( ) , , , 0 k b x y z t ρ = , ( ) 0 , , 0 k u x y z = and ( ) 0 , , 0 k v x y z = . Under these two assumptions, equation (3.1) becomes ( ) ( ) ( ) ( ) ( ) * * , , , d , , , d , k kl kl k kl k p t u t p t u t u t κ Γ Γ ∗ Γ = ∗ Γ + ∫ ∫ x x ξ x ξ x ξ (3.6) where kl kl κ δ = if ∈Ω ξ and 0 kl κ = if ∉Ω ξ . The integral equation (3.6) relates the displacements of a point ξ of the domain Ω with the displacements and tractions at the boundary Γ . In the frequency domain, equation (3.6) becomes ( ) ( ) ( ) ( ) ( ) * * , , , d , , , d , k kl kl k kl k p u p u u ω ω ω ω κ ω Γ Γ Γ = Γ + ∫ ∫ x x ξ x ξ x ξ ɶ ɶ ɶ ɶ ɶ (3.7) where a tilde over the variables denotes their Fourier transform with respect to time. Hence, ( ) * , , kl u ω x ξ ɶ and ( ) * , , kl p ω x ξ ɶ represent the three-dimensional fundamental displacements and tractions in the frequency domain, which can be calculated using the Thin Layer Method (TLM) (Kausel, 1981). 3.2.3 Regularization of the integral equation The fundamental displacements and fundamental traction are singular at the collocation point ξ and consequently equations (3.6) and (3.7) are not valid if ξ belongs to the boundary Γ . Two approaches can be followed to deal with this problem: a limiting process in which a spherical portion of the domain with radius tending to zero is excluded (or included) around the collocation ξ (Figure 3.1) (Dominguez, 1993); and a regularization procedure in which the singularities of the fundamental solutions are removed (François, 2008). Chapter 3 – Numerical tools for soil-structure interaction 66 Figure 3.3: Rigid footing resting on a half-space: boundary element mesh Next, to validate the procedure, the horizontal compliance ( HH x C u Ga = , considering a unit horizontal load) and the vertical compliance ( VV z C u Ga = , considering a unit vertical load) of a square footing ( a b = ) are compared with the compliances obtained by Wong and Luco (1978), for the dimensionless frequency 0 s a a C ω = ( s C G ρ = ) varying from 0 to 4. The size of the square foundation is 2 60 a = and the soil properties are 6 3.315 10 G= × , 4 2.82 10 ρ − = × and 1 3 ν = (consistent units). The soil-structure interface is divided into 100 equally sized squares, and within each square the displacement and traction fields are assumed to be constant (constant boundary elements). The fundamental solutions of the halfspace are obtained via the TLM with a model consisting of a elastic layer with thickness 5 0.1 s s h C π ω λ = = (divided into 25 thin layers of quadratic expansion) that rests on a PML with parameters 2 m = , 3 η = , 2 n = , 15 N = and 12 Ω = (see Chapter 2, equation 2.108). The elastic layer is added to the TLM model to improve the quality of the fundamental solutions near the source (recall that the TLM is a discrete method and so, near a point source, where the displacement and stress fields vary rapidly, thinner meshes are required in order to obtain a good accuracy). Figures 3.4 and 3.5 plot the results herein obtained and the results reported by Wong and Luco. Though not perfect, the agreement is good, which validates the procedure. Figure 3.4: Horizontal compliance HH C x y z 2 a 2 b 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.05 0.1 0.15 0.2 C HH a 0 Re(C HH ) Im(-C HH ) Re(C HH ) : Wong and Luco (1978) Im(-C HH ) : Wong and Luco (1978) Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 67 Figure 3.5: Vertical compliance VV C Note: In this example, because the boundary elements are constant, square shaped, and resting at the surface of the half-space, the calculation of the boundary element matrices ( ) ( ) I II ,I II U of equation (3.22) can be accomplished with a simplified procedure: instead of integrating the displacements induced by a point load at the i th boundary node on the surface of the j th boundary element, the 3x3 sub-matrix associated with the i th boundary node (rows) and j th boundary element (columns) can be taken as the displacements (all nine components) at the i th boundary node induced by a disk load applied at the j th boundary element. The radius of the disk load must be such that the area of the disk is the same as the area of the boundary element. This simplified procedure yields very good approximations as long as the boundary elements are horizontal, fairly square and of constant expansion. If the elements are not at the surface, this procedure cannot be used because it is not valid for the traction integrals. The displacements induced by disk loads are given in the work (Kausel and Peek, 1982b) and transcribed in Appendix IV. 3.2.6 Weak coupling – response to incoming wave fields As mentioned in the beginning of this chapter, in this work the propagation and the reception stages are decoupled, i.e., the wave-field that the track transmits to the soil is calculated disregarding the presence of the building in the far field. Hence, there is a weak coupling between the two sub-structures: the connection between the track and the structure is kept, while the connection between the structure and the track is not (Figure 3.6). So, consider that the tractions t p at the boundary t Γ (that result from the interaction between the track and the soil, and that are calculated without accounting for the removal of the volume of soil s Ω ) induce at the boundary s Γ the displacement field 0 u and the traction field 0 p (Figure 3.6). The displacements 0 u are calculated by placing collocation points on the boundary s Γ and then using equation (3.7) with t Γ = Γ . The traction field 0 p is obtained using the derivatives of 0 u in the strain-stress relations (2.6). The derivatives of 0 u are obtained by deriving equation (3.7) with respect to ξ — see (Dominguez, 1993) for details. 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.05 0.1 C VV a 0 Re(C VV ) Im(-C VV ) Re(C VV ) : Wong and Luco (1978) Im(-C VV ) : Wong and Luco (1978) Chapter 3 – Numerical tools for soil-structure interaction 68 Figure 3.6: Weak coupling between track and structure The first step to obtain the response of s Γ is to make it traction free, a condition that is needed to account for the volume s Ω of soil that is excavated. Such is accomplished by applying the tractions 0 − p at s Γ (same magnitude but opposite direction). These tractions induce an extra displacement field p u at s Γ calculated with ( ) 1 0 s s − = − + p u P I U p (3.35) where s P and s U are the boundary element matrices of the discretized boundary s Γ . The displacement field inc u that t p induces at s Γ (and that accounts for the volume of excavated soil s Ω ) is then ( ) 1 inc 0 0 0 s s − = + = − + p u u u u P I U p (3.36) (for s Γ at the surface, 0 p is null and so inc 0 = u u ). The second step to obtain the response of s Γ is to establish the compatibility of displacements and equilibrium of forces between the soil and the structure. The incident displacement field inc u needs to be taken into account in the compatibility equation (3.24), which becomes F B inc I I I = + N u u u ɶ ɶ ɶ (3.37) The equilibrium equation does not change. After combining the compatibility equation (3.37) and the equilibrium equation (3.25), a system of equations with the same form as equation (3.31) is obtained, with BEM F being calculated with t p s p 0 0 , u p t Γ s Γ s Ω Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 69 1 1 inc BEM II − − = − F TA C p TA B u ɶ (3.38) and with the remaining matrices and vectors being calculated as explained in section 3.2.5 for s Γ = Γ . Example: two weakly coupled rigid footings resting on a half-space Two square rigid footings ( 2 2 60 a b = = ) are placed at a distance 4 120 d a = = from each other (center to center). Both footings have concentrated masses placed at their center of gravity ( 500 M = ), which results in inertial forces in the translational degrees of freedom. The footings rest on a half-space whose properties are the same of those of the previous example. One of the foundations (footing 1) is loaded vertically and the vertical response of both footings ( VV,1 ,1 z C u Ga = and VV,2 ,2 z C u Ga = ) is calculated considering full coupling and weak coupling between the two foundations. Figures 3.7 and 3.8 compare the responses obtained considering the two coupling schemes for the dimensionless frequency 0 a varying between 0 and 1. The presence of the concentrated masses in the footings significantly modifies the behavior of the footing-soil system, as can be concluded from the comparison between Figures 3.5 and 3.7: both the shape and the magnitude of the response are different. It is also observed that the two coupling schemes yield different results. Nevertheless, the response of the weak coupling scheme follows the main trends of the response of the full coupling scheme. Figure 3.7: Vertical displacement of loaded footing: full coupling versus weak coupling 0 0.2 0.4 0.6 0.8 1 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 a 0 C VV,1 Re(C VV,1 ) - Full coupling Im(-C VV,1 ) - Full coupling Re(C VV,1 ) - Weak coupling Im(-C VV,1 ) - Weak coupli ng Chapter 3 – Numerical tools for soil-structure interaction 70 Figure 3.8: Vertical displacement of free footing: full coupling versus weak coupling 3.2.7 Final considerations In this section of chapter 3, the 3D BEM was introduced and coupled to the 3D FEM in order to perform dynamic analyses of structures interacting with the soil. To validate the implemented methodology, the example of a rigid footing resting on a half-space was solved and the results were compared with results available in the literature. Then, the definition of weak coupling between two structures was introduced and an example was shown where the results obtained considering weak coupling and full coupling were compared. It was observed that the two coupling approaches yielded different results, but that the major trends of the responses were kept. In this work, the 3D BEM-FEM methodology is used in the context of railway induced vibrations to analyze the response of structures to incoming wave-fields. The wave-fields are calculated using a 2.5D BEM-FEM procedure, which is explained in the following subsections. In the 2.5D BEM-FEM procedure, the presence of the building in the far field is disregarded, i.e., it is considered that the track and the buildings are weakly coupled. The consideration of buried structures is not attempted in this work, but some comments are proffered next about this issue: 1. The calculation of the boundary matrices for non-horizontal boundaries using the TLM becomes much more complicated than to integrate the fundamental displacements on a horizontal surface or to use the simplified procedure explained at the end of section 3.2.5. The difficulties arise because the boundary integrals cannot be calculated directly (unlike in the 2.5D case, explained in a later section of this chapter), which leads to the necessity of using special integration schemes (Gaussian integration, for example) and consequently an increase in the time needed for the computation of these integrals. Nevertheless, in the work (Kausel and Peek, 1982a) the BEM is formulated both in 2D and 3D spaces using the TLM, and these formulation can be used for the analyses of buried structures. 2. Alternative approaches can also be used: in recent years, some authors suggested the use of the Perfectly Matched Layer together with finite elements to calculate 0 0.2 0.4 0.6 0.8 1 -0.6 -0.4 -0.2 0 0.2 0.4 a 0 C VV,2 Re(C VV,2 ) - Full coupling Im(-C VV,2 ) - Full coupling Re(C VV,2 ) - Weak coupling Im(-C VV,2 ) - Weak coupling Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 71 the stiffness matrix s BEM ( ) Γ = K K of the soil interface (Basu and Chopra, 2003; Harari and Albocher, 2006). However, there are reports of incorrect results of the FEM-PML when the domain is layered, and that occurs due to grazing incidence of waves (Komatitsch and Martin, 2007). In this work, during some tests with FEM-PML to simulate a single stratum, for certain frequencies, the response of the stratum diverged from the expected. 3. Another alternative for the calculation of the stiffness matrix of the soil is the “SASSI” approach (Lysmer et al., 1999). In this approach, a flexibility matrix Ω F that relates forces and displacements of a grid of nodes that delimit the volume of soil to be excavated is calculated using the TLM. The flexibility matrix is then inverted, being thus obtained a stiffness matrix ( 1 − Ω Ω = K F ). The stiffness matrix s Γ K of the soil is finally obtained by subtracting from Ω K the stiffness s Ω K of the volume to be excavated, which is calculated with the FEM ( s s Γ Ω Ω = − K K K ). 4. As a final comment, independently of the approach that is followed, there is always a huge computational cost that cannot be avoided because of the need for a 3D mesh. 3.3 2.5D Finite Element Method 3.3.1 Introduction When the geometry of the structure is invariant in one of the directions, space-wavenumber- frequency domain (2.5D) analyses are usually more advantageous than space domain (3D) analyses. In essence, a 2.5D analysis consists in performing a Fourier transform of the longitudinal coordinate, which results in the reduction of the dimensionality of the problem by one, and consequently in the reduction of the 3D analysis to the summation of a set of 2D analyses. In the context of railway induced vibrations, since in most cases it can be assumed that the geometry of the track and the profile of the soil are invariant in the longitudinal direction, 2.5D analyses can be used to calculate the wave-fields that propagate in the track and soil and reach the buildings. The 2.5D analyses have been used by several authors in this context, as mentioned in the literature review presented in chapter 1. In this work, the track is modeled using 2.5D finite elements while the soil is modeled using 2.5D boundary elements. The two procedures are coupled in order to solve the track-soil system and to calculate the vibration field that propagates in the track and soil. Coupled 2.5D BEM-FEM schemes can be found in the literature, for example in Galvín et al. (2010) and in Sheng et al. (2005). The first work uses layered half-space fundamental solutions calculated with the stiffness matrices while the latter uses the analytical full-space fundamental solutions that are given in the same reference. The fundamental solutions used in the present work are the 2.5D fundamental solutions explained in chapter 2. 3.3.2 2.5D Finite Element Method Consider an elastic body that extends to infinity in one direction (longitudinal direction y ) and whose cross section S Ω is constant (Figure 3.9). It is assumed that the material properties are invariant in the y direction while within the cross section the material properties may vary in a stepwise fashion. Chapter 3 – Numerical tools for soil-structure interaction 72 Figure 3.9: Elastic solid: definition of invariance and angles θ Under the mentioned conditions, and following the same steps as for the TLM — see chapter 2, equations (2.2) to (2.13) — the wave equation in the elastic body can be written as (the following variables have the same meaning as in chapter 2) ( ) ( ) ( ) 2 2 2 2 2 2 2 2 2 xx xy yx xz zx yy yz zy zz x x y x z y y z z ρ  ∂ ∂ ∂ − + + + +  ∂ ∂ ∂ ∂ ∂   ∂ ∂ ∂ + + + + =  ∂ ∂ ∂ ∂  u D D D D D D D D D u b ɺɺ (3.39) On the other hand, the internal stresses s at any plane parallel to y can be calculated with ( ) T T cos sin x z θ θ = + s L L DLu (3.40) where θ is the angle of the outwards normal to the plane with respect to the x axis. At the interface between different materials, the internal stresses must be in equilibrium, i.e., 1 2 1 2 , θ π θ = = + s s (3.41) while at the boundary of the cross section S Ω the internal stresses must balance the external tractions t , i.e., = s t (3.42) As in the TLM case, one proceeds to discretize the domain to solve the wave equation. However, instead of discretizing only in the vertical direction z , in this case the domain is discretized in the xz plane. Hence, after dividing the cross section S Ω into plane finite elements, the displacement field in the solid is approximated by ( ) ( ) ( ) , , , x y z x z y =u N U (3.43) where N is a matrix containing the interpolation functions and U is a vector containing the displacements of the associated nodes. The interpolation functions N used in the 2.5D case are the same functions that are commonly used in plane strain or plane stress problems. After inserting the approximation (3.43) in the wave equation (3.39) and boundary conditions (3.41) and (3.42), it can be verified that these equations are not rigorously satisfied, due to the x y z 1 θ 2 θ 1 n 2 n 1 2 1 1 1 , , G ρ ν 2 2 2 , , G ρ ν Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 73 presence of unbalanced body forces and tractions. The discrete wave equation is derived after the application of the method of the weighted residuals and by requiring the virtual work done by the unbalanced forces within each finite element to be null, which results in the single finite element equation 2 2 xx xy xz yy zy zz y y y ∂ ∂ ∂ = + − + − − + ∂ ∂ ∂ U U U F MU G U B G U A B G U ɺɺ (3.44) The vector F contains the consistent external forces at the nodes of the finite elements (which result from the external tractions t and the body loads b ), while the finite element matrices M , yy A , y α B , αα G and xz G are given by ( , x z α = ) S T d d x z ρ Ω = ∫∫ M N N (3.45) S T d d yy yy x z Ω = ∫∫ A N D N (3.46) S S T d d d d T y y y x z x z α α α α α Ω Ω = − + ∫∫ ∫∫ B N D N N D N (3.47) S T d d x z αα α αα α Ω = ∫∫ G N D N (3.48) S S T T d d d d xz x xz z z zx x x z x z Ω Ω = + ∫∫ ∫∫ G N D N N D N (3.49) being xx ∂ ∂ = N N and zz ∂ ∂ = N N . After the assembly of the finite element matrices, one is left with a global system of differential equations with the same form as (3.44). To solve it, the displacements U and forces F are transformed to the wavenumber-frequency domain by means of the double Fourier transformations ( ) ( ) ( ) i , , e d d y t k y y k y t y t ω ω +∞ +∞ − − −∞ −∞ = ∫ ∫ U U (3.50) ( ) ( ) ( ) i , , e d d y t k y y k y t y t ω ω +∞ +∞ − − −∞ −∞ = ∫ ∫ F F (3.51) In the transformed domain, the system (3.44) changes into ( ) ( ) 2 2 i yy y y xy zy xx xz zz k k ω   = + + + + + −   F A B B G G G M U (3.52) All matrices in (3.52) are symmetric, except for xy B and zy B , which are skew-symmetric. However, for cross-anisotropic materials whose constitutive matrix is given in equation (2.1) of chapter 2, it is possible to transform these matrices into symmetric matrices by means of a similarity transformation that consists in multiplying the rows of (3.52) that are related to the y direction by i − and the columns related to the y direction by i . The referred to rows and columns have indexes 2 3 l j = + , with 0,1,2,... j = . This transformation solely affects the matrices xy B and zy B and the vectors F and U . After the transformation, the system (3.52) becomes Chapter 3 – Numerical tools for soil-structure interaction 74 ( ) ( ) 2 2 yy y y xy zy xx xz zz k k ω   = + + + + + −   f A B B G G G M u (3.53) where f and u are obtained from F and U after multiplying the rows l by i − . Also, xy B and zy B are obtained from xy B and zy B by reversing the sign of the columns with index l . After solving the system of equations (3.53) for u , the displacements U in the wavenumberfrequency domain are recovered by multiplying every row l of u by +i . [Note1: the above description considered the cross-anisotropic material defined in eq. (2.1). In some cases, it is more convenient to define materials whose direction of anisotropy is the longitudinal direction instead of the vertical one. For example, it is common to model the sleepers (which are discontinuous) as a continuous anisotropic slab, where the in-plane behavior differs from the out-of-plane behavior (Alves Costa, 2011). In that case, the constitutive matrix ( ) ( ) ( ) 1 1 1 1 2 1 2 1 2 1 000 000 000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 xyz xz xz xz xz xyz xyz xz y xz xyz xz xz xz xz y y xz xz y y E E E E E E E E E E E E νν ν ν ν ν ν ν ν − + + +   − −     − −     − −   =               D (3.54) is more appropriate. xz E and y E are the in-plane and out-of-plane elastic modulus, xz ν and y ν are the in-plane and out-of-plane Poisson’s ratio, and xyz ν is the cross Poisson’s ratio. This last variable can assume negative values and must be such that matrix D is definite positive. The equations presented in this section are still valid for this constitutive matrix.] [Note2: the above formulation considers only volume elements. In order to consider different elements, such as beams or shells, different differential equations and discretizations have to be used, resulting in a final system of equations that will contain terms in 3 y k and 4 y k (Galvín et al., 2010). For an Euler beam, the system of equations is ( ) ( ) ( ) 4 2 2 2 4 2 x x y x y y y z z y z f EI k A u f EAk A u f EI k A u ω ρ ω ρ ω ρ = −  = −  = −   (3.55) where E is the Young’s modulus of the beam, ρ is the mass density, , x z I are the moments of intertia, A is the cross section area, , , x y z f represent the external forces and , , x y z u represent the displacements of the axis of the beam.] 3.3.3 Example - dispersion curves of a UCI861-3 rail To validate the implementation of the procedure presented in this section, the dispersion curves of a UCI861-3 rail are calculated and compared with the curves obtained by Gavric (1995), who also used a 2.5D FEM procedure. Dispersion curves correspond to the pairs Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 75 ( , y k ω ) that result in singular systems of equations (3.52) or (3.53). In other words, the dispersion curves plot the values y k that are solutions of ( ) ( ) 2 2 det 0 yy y y xy zy xx xz zz k k ω   + + + + + − =   A B B G G G M (3.56) as function of the frequency ω . The solutions of (3.56) may be real or complex: the real solutions correspond to waves that propagate while the complex solutions correspond to waves that evanesce (i.e., attenuate) with the distance to the source. The 2.5D FEM methodology presented in this section is employed to determine the dispersion curves of the propagating modes of the rail. With that intention, the rail section is divided into 166 quadrilateral elements and a total of 214 nodes, resulting in the mesh shown in Figure 3.10. For that mesh, the global matrices yy A , y α B , αβ G and M are calculated and then used in equation (3.56). Figure 3.10: Used mesh for the rail section (dimensions in mm) In Figure 3.11, the results obtained in this work (black) are compared with the curves obtained by Gavric (red). The agreement of the curves is good namely for the lower frequencies. For the higher frequencies (above 3000 Hz), even though the shapes of the curves are the same, there is a shift between the results obtained in this work and the results obtained by Gavric. The disagreement might be justified by discrepancies in the material parameters: while in this work the material properties are Young’s Modulus 200 GPa E = , mass density 3 7859 kg / m ρ = , and Poisson’s ratio 0.28 ν = , in (Gavric, 1995) the material properties are not specified. [Note: equation (3.56) is solved for y k by determining the eigenvalues and eigenvectors of ( ) ( ) 2 2 yy y y xy zy xx xz zz k k ω   + + + + + − =   A B B G G G M 0 f (3.57) If the rows 2 3 l j = + ( 0,1,2,... j = ) of (3.57) are multiplied by y k and the columns l are divided by y k , the symmetric eigenvalue problem (3.57) can be reduced to a non-symmetric general eigenvalue problem of the form ( ) 2 y k + = A C 0 f (3.58) Chapter 3 – Numerical tools for soil-structure interaction 82 ( ) ( 1) ( 1) ( 2) ( 2) i2 2 2 2 2 p p p nj nj nj p p nj nj l l l J I x I x l l I x I x ξ ξ ξ ξ − − − −       = − − + − + +                   − − − − +             (3.84) The integrals ( ) ( 2) nj I x − required to evaluate the coefficients ( ) ej kl U and ( ) ej kl P are given in Table 3.2. Case 3: ( ) ( ) 2 2 e j S x S x x = = The Fourier transform of ( ) 2 S x , according to (3.74), is ( ) 22i i i i i i i 22 2 2 2 2 2 22 3 2 i 2i e d e e e +e e e 4 x x x x x x x l k l k l k l k l k l k l k x x x x x l l l S k x x k k k − − − − −       = = − − + + −             ∫ ɶ (3.85) The coefficients ( ) ej kl U are then obtained with ( ) 2i i i i ( ) ( ) 2 2 2 2 2 i i 2 2 3 1 i e e e +e 2 4 2i e e d x x x x x x l l l l k x k x k x k x ej mn kl kl x x x l l k x k x x x l l U U k k k k k ξ ξ ξ ξ ξ ξ π         +∞ + − + −                 −∞     + −              =− − + +                  −         ∫ (3.86) and the coefficients ( ) ej kl P with ( ) 2i i i i ( ) 2 2 2 2 2 i i 2 2 3 1 i e e e +e 2 4 2i e e d x x x x x x l l l l k x k x k x k x ej mn kl kl x x x l l k x k x x x l l P t k k k k k ξ ξ ξ ξ ξ ξ π         +∞ + − + −                 −∞     + −              =− − + +                  −         ∫ (3.87) From equations (3.86) and (3.87) it can be concluded that the calculation of the coefficients ( ) ej kl U and ( ) ej kl P follow exactly the same steps as the calculation of ( ) ( )mn kl u x (section 2.5.1) and ( ) ( )mn kl t x (section 2.5.3), changing only the integrals ( ) ( ) p nj I x by the integrals ( ) p nj J of the form 2 ( ) ( 1) ( 1) ( 2) ( 2) ( 3) ( 3) i4 2 2 2 2 2i 2 2 p p p nj nj nj p p nj nj p p nj nj l l l J I x I x l l l I x I x l l I x I x ξ ξ ξ ξ ξ ξ − − − − − −       = − − − − − + +                   − − + − + +                   − − − − +             (3.88) The integrals ( ) ( 3) nj I x − required to evaluate the coefficients ( ) ej kl U and ( ) ej kl P are given in Table 3.3. Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 83 Table 3.1: Closed form expressions for ( 1) nj I − ( 2 2 Im 0 j y k k − < ) ( ) ( ) ( ) ( ) 2 2 2 2 i i ( 1) 1 1 1 2 2 i i ( 1) 1 2 2 22 2 i i ( 1) 1 3 3 2 1 2 1 2 1 2 sign e d 1 e 2i 1 e d sign e i e 2 sign e d e e 2i j y x y j y x y x k k x k x j x j x y j k x k k x y k x j x j x y jj y k x k k x j x j x j x I k K k k k k I k K k k kk k x I k K k k π π π +∞ − − − − − −∞ +∞ − − − − − − −∞ +∞ − − − − − −∞   = = −     −     = = − +   −   = = − − ∫ ∫ ∫ ( ) ( ) ( ) 2 2 2 2 2 2 2i i ( 1) 1 4 4 2 2 2 2 2 2 i i ( 1) 1 5 5 2 2 i ( 1) 1 6 6 1 2 1 2 1 2 sign 1 e e d e 2i 1 e d e 2 i sign e j y y j y x j y x x k x k x k k x y k x j x j x y j j j y j k k x k x j x j x j j y k x j x j x k x I k K k k k k k k k I k K k k k k x k I k K dk π π π − − +∞ − − − − − −∞ +∞ − − − − − −∞ +∞ − − − −∞           = = + −   −−     = = − = = ∫ ∫ ∫ ( ) 2 2 i 2 2 1 e 2i j y k k x y j y j k k k − −   −     − Table 3.2: Closed form expressions for ( 2) nj I − ( 2 2 Im 0 j y k k − < ) ( ) ( ) ( ) ( ) 2 2 2 2 i i ( 2) 2 1 1 2 2 2 2 i i ( 2) 2 2 2 2 2 2 2 2 2 2 2 i ( 2) 2 3 3 1 2 1 2 1 2 1 e e d i 2 sign 1 e e e d 2i e d j y x y j y x x k k x k x j x j x y j j y k x k k x y k x j x j x y j y y j j y j k x j x j x I k K k x k k k k k x I k K k k k k k k k k k I k K k π π π − − +∞ − − − −∞ − − − +∞ − − − −∞ +∞ − − − −∞     = = − −   −−       = = + −   − −     = ∫ ∫ ( ) ( ) ( ) 2 2 2 2 2 i 22 2 i 2 i ( 2) 2 4 4 2 2 22 2 2 2 2 i i ( 2) 2 5 5 2 2 1 2 1 2 1 e e i 2 e 1 e e d i 2 sign e d 1 e 2i y j y j y y x j x k x k k x jyj y k k x k x y k x j x j x y j y j j y j j y k k k x j x j x j y j kkk k k x I k K k k k k k k k k k k x I k K k k k k π π − − − − − − +∞ − − − −∞ +∞ − − − − − −∞     = − +   −     −   = = + +   −− −     = = − − ∫ ∫ ∫ ( ) 2 2 2 i i 2 2 6 6 2 2 2 2 1 2 e e d i 2 y j y x x k k x y k x j x j x j y j j y k I k K k x k k k k k π − − +∞ − − − −∞           = = − −   −−   ∫ Chapter 3 – Numerical tools for soil-structure interaction 84 Table 3.3: Closed form expressions for ( 3) nj I − ( 2 2 Im 0 j y k k − < ) ( ) ( ) ( ) 2 2 2 2 i 2 i ( 3) 3 1 1 2 2 2 2 2 2 i 2 i ( 3) 3 2 2 2 2 22 2 2 2 2 ( 3) 3 1 2 1 2 1 2 sign 1 e e d 2 2i e 1 e e d i 2 j y x j y y x k k x k x j x j x y j y j y j k k x k x y k x j x j x y y j y j j y j j y j x xx I k K k k k k k k k k x I k K k k k k k k k k k k k I k π π π − − +∞ − − − −∞ − − − +∞ − − − −∞ −     = = − + −  − − −    −   = = + +   −− −     = ∫ ∫ ( ) ( ) ( ) ( ) ( ) ( ) 2 2 2 2 i 2 i 3 32 2 2 2 2 2 2 i 2 2 2 2 i ( 3) 3 4 4 2 2 2 2 2 2 2 2 2 2 2 2 1 2 e sign 1 e e d 2i 2 e sign e e d 2i 2 j y y x j y y x k k x k x y k x j x y y j j j y j k k x k x y j y k x j x j x y j y j y y j j y j k x K k k k k k k k k k k k xx I k K k k k k k k k k k k k π − − − +∞ − − −∞ − − − +∞ − − − −∞     = + −   −−       −   = = − + + −   −− −    ∫ ∫ ( ) ( ) ( ) 2 2 2 2 i i ( 3) 3 5 5 2 2 2 2 i 2 i ( 3) 3 6 6 2 2 2 2 2 2 1 2 1 2 1 e e d i 2 sign 1 e e d 2 2i j y x j y x k k x k x j x j x j y j j y k k x y k x j x j x y j y j j y j I k K k x k k k k k x k x I k K k k k k k k k k π π − − +∞ − − − −∞ − − +∞ − − − −∞      = = − −   −−       = = − + −   − − −  ∫ ∫ Considerations concerning horizontal boundaries The calculation of the coefficients ( ) ej kl U involves only the components of the modal shapes at the elevation of the collocation point and at the elevation of the boundary element. By contrast, the calculation of the coefficients ( ) ej kl P involves the components of all TLM nodes that compose the thin-layer delimiting the boundary element. Since the boundary elements are placed at the interface between two consecutive thin-layers, a decision is required as to whether to consider the upper or the lower thin-layer. When the collocation point is not contained in the boundary element, it is immaterial which thin-layer is used. On the other hand, when the collocation point is contained in the boundary element, the value of ( ) ej kl P depends on the thin-layer selected for the evaluation. The rule used in this work is that if the outwards normal faces up, the thin-layer located below the boundary is employed in the calculation of ( ) ej kl P , otherwise the thin-layer above is used. By following this procedure, collocation points on horizontal boundaries are circumvented as depicted in Figure 3.13. It is important to note that according to this procedure, when the collocation point ξ is at an edge of a boundary element ( 2 x l ξ = ± ), then the coefficient ( ) ej kl P is calculated considering that the boundary is distorted as shown in Figure 3.14. This aspect is important for the treatment of corners, i.e., points where horizontal boundaries meet vertical boundaries (Figure 3.13b-c). Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 85 Figure 3.13: Exclusions in the domain at the collocation points. a) smooth horizontal boundary; b) concave corner; c) convex corner Figure 3.14: Deflection of the horizontal boundary when the collocation point is at one extreme Due to the homogeneity of layered domains in the horizontal direction, the fundamental solutions depend on the horizontal distance between the source and the receiver, and not on their absolute horizontal coordinates. Hence, if the horizontal boundaries are dicretized in such a way that the boundary elements have the same length and expansion order, then the boundary coefficients ( ) ej kl U and ( ) ej kl P can be reused whenever the distance between the collocation points and the boundary elements is repeated. This fact can reduce significantly the cost of computation of matrices U and P . Validation Consider a homogeneous full-space with mass density 1 ρ = , shear modulus 1 G = , Poisson’s ratio 0.25 ν = and hysteretic damping 0.001 p s ξ ξ = = . The full-space is simulated with a TLM model consisting of an elastic layer with thickness 2 H = that is divided into 40 thin-layers of quadratic expansion and that is supplemented with two PMLs, one at the top and the other at the bottom, with parameters 2 m = , 2 η = , 8 Ω = and 10 N = (see section 2.6). Consider also a horizontal boundary element (of quadratic expansion) whose width is 0.1 BEM l= and place it at the depth 0 BEM z = and horizontal coordinate 0 BEM x = . The shape functions associated with the nodes of the boundary elements are (from left to right) n ξ upper thin-layer lower thin-layer boundary deflected boundary ε a) b) c) Chapter 3 – Numerical tools for soil-structure interaction 86 ( ) ( ) ( ) 2 12 2 22 2 32 2 4 1 2 BEM BEM BEM BEM BEM x x S x l l x S x l x x S x l l = − +    = −    = +   (3.89) In the following, the coefficients kl U and kl P associated with the middle node (node 2) are calculated with the procedure described in this section and compared with the results obtained through simple numerical integration (1000 integrating points) of the analytical solutions (Tadeu and Kausel, 2000). Three collocation points are considered: the first is placed at the position ( ) 1 0.5, 0.5 = ξ (outside the boundary element), the second at ( ) 2 0.05, 0 = ξ (right edge of the boundary element) and the third at ( ) 3 0, 0 = ξ (center of the boundary element). For the last collocation point and for the integration of the analytical solutions, the boundary contour is deflected as indicated in Figure 3.13a in order to avoid the singularities. The assumed frequency is 1Hz f = ( 2 rad/s ω π = ). Figure 3.15-3.17 plot the comparison between the results obtained with the two approaches. As can be observed from Figure 3.15, the two approaches yield practically the same results, thus validating the procedure. Figure 3.16 provides exactly the same conclusions: discrepancies can be noticed for some components of kl U and kl P but these differences exist merely due to some residual values of the TLM results. In Figure 3.17, the agreement is also very good: also in this example, the differences are caused by residual values. Notice that the last two collocation points are placed inside the boundary element and that no special treatment is given to the fundamental solutions of the TLM and that the results are still correct. This confirms that using this procedure the coefficients kl c are already accounted for during the calculation of the boundary integrals kl P : for the collocation point 3 ξ , the components xx P , yy P and zz P correspond to 0.5 kl kl c δ = (smooth boundary). Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 87 Figure 3.15: Boundary coefficients kl U and kl P for 1 ξ . Solid lines – TLM solution; circles – analytical solutions. Blue – real component; red – imaginary component 0 10 20 -0.02 0 0.02 k y U xx 0 10 20 -5 0 5x 10 -3 k y U xy 0 10 20 -2 0 2 4x 10 - 3 k y U xz 0 10 20 -5 0 5x 10 -3 k y U yx 0 10 20 -5 0 5 10 x 10 - 3 k y U yy 0 10 20 -5 0 5x 10 - 3 k y U yz 0 10 20 -2 0 2 4x 10 -3 k y U zx 0 10 20 -5 0 5x 10 -3 k y U zy 0 10 20 -0.02 0 0.02 k y U zz 0 10 20 -0.02 0 0.02 k y P xx 0 10 20 -0.02 0 0.02 k y P xy 0 10 20 -0.02 0 0.02 k y P xz 0 10 20 -0.02 0 0.02 k y P yx 0 10 20 -0.02 0 0.02 0.04 k y P yy 0 10 20 -0.2 0 0.2 k y P yz 0 10 20 -0.04 -0.02 0 0.02 k y P zx 0 10 20 -0.04 -0.02 0 0.02 k y P zy 0 10 20 -0.05 0 0.05 k y P zz Chapter 3 – Numerical tools for soil-structure interaction 88 Figure 3.16: Boundary coefficients kl U and kl P for 2 ξ . Solid lines – TLM solution; circles – analytical solutions. Blue – real component; red – imaginary component 0 10 20 -0.05 0 0.05 k y U xx 0 10 20 -3 -2 -1 0x 10 -3 k y U xy 0 10 20 -1 0 1x 10 -13 k y U xz 0 10 20 -3 -2 -1 0x 10 -3 k y U yx 0 10 20 -0.02 0 0.02 0.04 k y U yy 0 10 20 -5 0 5x 10 -14 k y U yz 0 10 20 -1 0 1x 10 -13 k y U zx 0 10 20 -5 0 5x 10 -14 k y U zy 0 10 20 -0.05 0 0.05 k y U zz 0 10 20 -10 -5 0 5x 10 -12 k y P xx 0 10 20 -2 0 2 4x 10 -12 k y P xy 0 10 20 -0.2 0 0.2 k y P xz 0 10 20 -2 0 2 4x 10 - 12 k y P yx 0 10 20 -10 -5 0 5x 10 -12 k y P yy 0 10 20 -0.4 -0.2 0 k y P yz 0 10 20 -0.2 0 0.2 k y P zx 0 10 20 -0.1 0 0.1 k y P zy 0 10 20 -5 0 5x 10 -12 k y P zz Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 89 Figure 3.17: Boundary coefficients kl U and kl P for 3 ξ . Solid lines – TLM solution; circles – analytical solutions. Blue – real component; red – imaginary component 0 10 20 - 0.1 0 0.1 k y U xx 0 10 20 -2 0 2 x 10 -19 ky U xy 0 10 20 -1 0 1 ky U xz 0 10 20 -2 0 2 x 10 -19 k y U yx 0 10 20 -0.02 0 0.02 0.04 k y U yy 0 10 20 -1 - 0.5 0x 10 -13 k y Uyz 0 10 20 -1 0 1 k y U zx 0 10 20 0 0.5 1 x 10 -13 k y U zy 0 10 20 -0.05 0 0.05 k y U zz 0 10 20 -0.5 0 0.5 1 k y P xx 0 10 20 -1 0 1 k y P xy 0 10 20 -5 0 5 10 x 10 -10 k y P xz 0 10 20 -1 0 1 k y P yx 0 10 20 -0.5 0 0.5 1 k y P yy 0 10 20 -0.4 -0.2 0 k y P yz 0 10 20 -5 0 5 10 x 10 -10 k y P zx 0 10 20 0 0.1 0.2 k y P zy 0 10 20 -0.5 0 0.5 1 k y P zz Chapter 3 – Numerical tools for soil-structure interaction 90 3.4.4 Vertical boundary elements Vertical boundaries are defined by a constant horizontal coordinate BEM x (Figure 3.18). If it is assumed that the load is applied at the depth n z ( th n interface of the TLM model) and that the boundary element is placed between depths 1 m z and 2 m z ( th 1 m and th 2 m interfaces of the TLM model), then the integrals in equations (3.68) and (3.69) can be replaced by integrals of the form (for convenience, the variables , y k ω are dropped) ( ) ( ) ( ) 2 1 ( ) ( ) ( ) ( ) ( ) d m ej mn j kl kl BEM m e m m U u x x N z S z z ξ = = − ∑∫ (3.90) ( ) ( ) ( ) 2 1 ( ) ( ) ( ) ( ) ( ) d m ej mn e kl kxl BEM m j m m P x x N z S z z ξ σ = = ± − ∑∫ (3.91) In these equations, the factors ( ) ( ) ( ) ( ) mn kl BEM m u x x N z ξ − and ( ) ( ) ( ) ( ) mn kxl BEM m x x N z ξ σ − represent the vertically interpolated displacements and tractions fields, with ( ) ( ) m N z being the TLM shape function associated with the th m interface. In equation (3.91), the positive sign must be used if the outwards normal is in the positive x direction, while the negative sign must be used otherwise. Figure 3.18: Vertical boundary element Since ( ) ( )mn kl BEM u x x ξ − and ( ) ( )mn kxl BEM x x ξ σ − are nodal values and therefore do not depend on the depth z , the expressions (3.90) and (3.91) can be replaced by ( ) ( ) ( ) 2 1 ( ) ( ) ( ) ( ) ( ) d m ej mn j kl kl BEM m e m m U u x x N z S z z ξ = = − ∑∫ (3.92) ( ) ( ) ( ) 2 1 ( ) ( ) ( ) ( ) ( ) d m ej mn e kl kxl BEM m j m m P x x N z S z z ξ σ = = ± − ∑∫ (3.93) Thus, only the integrals of the form ( ) ( ) ( ) ( ) ( ) d e m j N z S z z ∫ need to be evaluated. Since ( ) ( ) m N z and ( ) ( ) ( ) e j S z are both polynomial functions, these integrals can be evaluated in closed form. As final note, since the displacements are interpolated in the vertical direction using polynomial functions, the singular behavior of the fundamental solutions is not captured. Hence, when the collocation point lies within the vertical boundary element, in the calculation of ( ) ej kl P the term kl c is not accounted for. Nonetheless, since the boundary elements are vertically oriented and the fundamental solutions are symmetric with respect to vertical planes, the resulting value for the missing term is 0.5 kl kl c δ = . In this way, for nodes that belong to vertical boundary elements and that do not correspond to corners, the term 2 BEM l ξ e Γ 1 m n 2 m 2 BEM l Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 91 0.5 kl kl c δ = must be added to the diagonal of P associated with the node. When the node corresponds to a corner, two situations occur: 1. Concave corner (Figure 3.13b): in this case, because the horizontal boundary element already accounts for the quarter of circle of the deflected boundary (Figure 3.14), then the factor kl c must only account for the remaining semi-circle, and so 0.5 kl kl c δ = ; 2. Convex corner (Figure 3.13c): the horizontal boundary element already accounts for the quarter circle of the deflected boundary (Figure 3.14), and so the factor kl c is null. Validation Recall the homogeneous full-space used to validate the horizontal boundary elements, which is simulated with the same TLM model, and consider a vertical boundary element with width 0.1 BEM l= , centered at ( ) ( ) , 0,0 BEM BEM x z = and of quadratic expansion (the boundary element is contained in two distinct thin-layers). The shape functions of the boundary element are the same as in equation (3.89), with the argument being replaced by z . In the following, the boundary coefficients kl U and kl P associated with the middle node (node 2) are calculated with the procedure described in this section and compared with the results obtained through simple numerical integration (1000 integrating points) of the analytical fundamental solutions (Tadeu and Kausel, 2000). Three collocation points are considered: the first is placed at the position ( ) 1 0.5, 0.5 = ξ (outside the boundary element), the second at ( ) 2 0, 0.05 = ξ (upper edge of the boundary element) and the third at ( ) 3 0, 0 = ξ (center of the boundary element). The assumed frequency is again 1Hz f = ( 2 rad/s ω π = ). Figure 3.19-3.21 plot the comparison between the results obtained with the two approaches. Such as in the case of the horizontal boundary element, a very good agreement is obtained between the two procedures. There is however a shift in the real components of the boundary coefficients xz P and zx P for 2 = ξ ξ and this difference does not vanish with the refinement of the TLM model. Nevertheless, as it is seen by the examples described in the next section, despite these differences, the results obtained with the TLM-BEM exhibit a good quality. Chapter 3 – Numerical tools for soil-structure interaction 98 Figure 3.22: Square tunnel inside a horizontally layered domain Since the cross section of the tunnel is rigid, the displacements of the walls of the tunnel can be described as function of the translation and rotation of the tunnel, i.e., ( ) ( ) ( ) Tunnel Tunnel Tunnel Tunnel Tunnel Tunnel Tunnel Tunnel , 1 0 0 0 0 , 0 1 0 0 , 0 0 1 0 0 x y x z y x z y z u u u x z z u u x z z x u x z x θ θ θ                 = = − =             −             Nu N u (3.102) On the other hand, the pressures that the layered domain transmits to the walls of the tunnel induce at the center of the tunnel forces and moments that are calculated with ( ) ( ) ( ) T Tunnel T , , d , x x y z x y z y z p x z f f f m m m p x z p x z Γ       = = Γ         ∫ f N (3.103) where Γ represents the boundary of the tunnel. After discretizing the surface of the layered domain that is in contact with the tunnel into boundary elements and N boundary nodes, the nodal displacements j u at the boundary are obtained with ( ) ( ) 1 1 1 Tunnel U U , , N N N x z x z         = =             u N N u N u N ⋮ ⋮ (3.104) The forces Tunnel f are obtained from the boundary pressures j p through , , , u u u u G ρ ν ξ 1 1 1 1 , , , G ρ ν ξ 2 2 2 2 , , , G ρ ν ξ , , , l l l l G ρ ν ξ 1 H 2 H 2 H L x z y 2 H Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 99 ( ) ( ) ( ) ( ) 1 Tunnel T T P P 1 , , d , , d N N x z S x z x z S x z Γ Γ       = = Γ Γ           ∫ ∫ p f N N N N p ⋮ ⋯ (3.105) with ( ) , j S x z being the shape function associated to the th j boundary node. Replacing eqs. (3.104) and (3.105) in equation (3.67) yields ( ) ( ) 1 Tunnel Tunnel P U − + =N U P C N u f (3.106) and so the compliance matrix corresponds to the 6 by 6 matrix F obtained with ( ) ( ) 1 1 P U − − = +F N U P C N (3.107) In the subsequent examples, the components of the compliance matrix F are evaluated using the 2.5D BEM methodology explained earlier. The tunnel is given the cross section 1[m] H L = = and each edge of the tunnel is divided into 5 boundary elements of quadratic expansion (3 nodes per boundary element). The total number of nodes is then 40 N = . To validate the results, the compliance matrices are also calculated using a finite element model coupled with PMLs, which are obtained as explained in (Kausel and Barbosa, 2011). The excitation frequency is 2 [rad s] ω π = and the wavenumbers y k range from 0 to 6 [rad m] π (301 wavenumbers). Homogeneous full-space The material properties of the full-space are: mass density 3 1 2 1 kg m u l ρ ρ ρ ρ = = = = ; shear modulus 1 2 1Pa u l G G G G= = = = ; Poison’s ratio 1 2 0.25 u l ν ν ν ν = = = = ; hysteretic damping 1 2 0.01 u l ξ ξ ξ ξ = = = = . The TLM model consists of the 4 macro-layers identified in Figure 3.22, where the upper and the lower semi-infinite elements are modeled with PMLs (with parameters 2, 8, 10, 2 N m η = Ω = = = ; see chapter 2 for definition of variables), and the middle layer satisfy 1 2 2m H H= = and are divided into 40 thin-layers of quadratic expansion. Due to symmetry conditions, only the components xx f , yy f , zz f , x x f θ θ , y y f θ θ , z z f θ θ , z z x x f f θ θ = − and x x z z f f θ θ = − are non-zero. Also, due to the geometry of the problem, xx zz f f = , x x z z f f θ θ θ θ = and z x x z f f θ θ = . Hence, considering only the five compliance components xx f , yy f , x x f θ θ , y y f θ θ and z x f θ , it is possible to describe the entire system. In Figure 3.23, the 5 components of the compliance matrix obtained with the proposed procedure (solid lines) are compared with the results obtained with the FEM (black circles). Blue is used for the representation of the real part, while red is used for the imaginary part. Figure 3.23 shows that the two approaches yield virtually identical results, leading to the conclusion that both procedures are correct. It can also be observed that the in-plane components ( xx f , x x f θ θ and z x f θ ) present singularities at 2 y S S k k C ω π = = = . It should be noted that, because in this example the soil is a homogeneous, infinite space, it follows that the classical BEM that uses the fundamental solutions of a full-space has a clear advantage over the use of BEM-TLM. However, this problem of very simple geometry is used solely for validation purposes. In the next examples, the use of full-space fundamental solutions requires the discretization not only of the edges of the tunnel but also of the free- Chapter 3 – Numerical tools for soil-structure interaction 100 surfaces and of the interface between different layers, and now the TLM offers clear advantages inasmuch as these interfaces need not to be discretized. Figure 3.23: Tunnel compliances for the full-space case. Solid lines = 2.5D BEM (real part – blue; imaginary part – red). Black circles = FEM Homogeneous layer free in space The free layer consists of the two intermediate macro layers depicted in Figure 3.22 ( 1 2 2 m H H= = ). The material properties of the free layer are the same of the fullspace considered in the previous example. The TLM model is similar to the one used therein, but with the upper and lower PMLs excluded. Again, due to symmetry conditions, only the components xx f , yy f , zz f , x x f θ θ , y y f θ θ , z z f θ θ , z z x x f f θ θ = − and x x z z f f θ θ = − do not vanish. However, the identities xx zz f f = , x x z z f f θ θ θ θ = and z x x z f f θ θ = do not hold, and so a total of eight components of the compliance matrix are needed to describe the system. Figure 3.24 shows the eight components obtained with the proposed methodology and with the FEM. Once again, the results obtained with the two procedures match perfectly. 0123 - 0.4 - 0.2 0 0.2 0.4 0 1 2 3 -0.1 -0.05 0 0.05 0.1 0123 - 0.4 - 0.2 0 0.2 0.4 0.6 0 1 2 3 -0.2 -0.1 0 0.1 0.2 0.3 xx f 2 y k π 2 y k π 2 y k π 2 y k π yy f x x f θ θ y y f θ θ 012 3 -0.4 -0.3 -0.2 -0.1 0 0.1 z x f θ 2 y k π Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 101 Figure 3.24: Tunnel compliances for the free layer in space. Solid lines = 2.5D BEM (real part – blue; imaginary part – red). Black circles = FEM 0 1 2 3 -0.3 -0.2 -0.1 0 0.1 0.2 0 1 2 3 -0.15 -0.1 -0.05 0 0.05 0.1 0 1 2 3 -0.2 -0.1 0 0.1 0.2 0 1 2 3 -1 -0.5 0 0.5 1 0 1 2 3 -0.4 -0.2 0 0.2 0.4 0 1 2 3 -0.4 -0.2 0 0.2 0.4 0 1 2 3 -0.4 -0.3 -0.2 -0.1 0 0.1 0 1 2 3 -0.1 0 0.1 0.2 0.3 xx f xx f x x f θ θ z x f θ xx f x x f θ θ x x f θ θ x z f θ 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π Chapter 3 – Numerical tools for soil-structure interaction 102 Homogeneous half-space The material properties of the homogeneous half-space are the same as in the previous case. The TLM model differs from the model in the first example in that the upper PML is excluded. In this case, 12 distinct components are needed to define the compliance matrix, whose structure is 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 y z x x x x x x y y y y z z y z z z xx x x yy yz y yz zz z y z x x f f f f f f f f f f f f f f f f f f θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ         −   =   −       − −     F (3.108) The 12 components of the compliance matrix are plotted in Figure 3.25. A good agreement is once again reached. Layered half-space The case of a non-homogeneous half-space is considered next. The properties of the layers, based on Figure 3.22, are the following: 0, 0 u u G ρ = = (the upper half-space does not exist) 3 1 1 1 1 1 1.2kg m , 1.0Pa, 0.25, 0.01, 2m G H ρ ν ξ = = = = = 3 2 2 2 2 2 1.3kg m , 2.0Pa, 0.3, 0.01, 2m G H ρ ν ξ = = = = = 2 2 2 2 , , , l l l l G G ρ ρ ν ν ξ ξ = = = = Each of the physical layers are modeled with 40 thin-layers based on a quadratic expansion. The lower half-space is modeled with PMLs with the same parameters used in the first example ( 2, 8, 10, 2 N m η = Ω = = = ). As in the case of the half-space, the 12 components given by eq. (3.108) are needed to define the compliance matrix F . These compliance components, obtained with the proposed procedure and with the FEM, are plotted in Figure 3.26. Again, the agreement between the proposed method and the FEM is very good. It can be concluded from this example that the BEM based on the TLM fundamental solutions can correctly simulate horizontally layered domains without the need to discretize the interfaces between layers, as is necessary in the standard BEM. Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 103 Figure 3.25: Tunnel compliances for the half-space. Solid lines = 2.5D BEM (real part – blue; imaginary part – red). Black circles = FEM zz f 0 1 2 3 -0.4 -0.2 0 0.2 0.4 0 1 2 3 -0.4 -0.2 0 0.2 0.4 0.6 0 1 2 3 -0.01 0 0.005 0.01 0 1 2 3 -0.04 -0.02 0 0.02 0.04 -0.005 0 1 2 3 -0.1 -0.05 0 0.05 0.1 0 1 2 3 -0.2 -0.1 0 0.1 0.2 0.3 0 1 2 3 -0.4 -0.3 -0.2 -0.1 0 0.1 0 1 2 3 -0.1 0 0.1 0.2 0.3 0.4 yy f y y f θ θ z z f θ θ 0 1 2 3 -0.4 -0.2 0 0.2 0.4 0 1 2 3 -0.4 -0.2 0 0.2 0.4 0.6 0 1 2 3 -0.02 -0.01 0 0.01 0.02 0 1 2 3 -0.04 -0.02 0 0.02 0.04 0.06 xx f x x f θ θ 2 y k π y x f θ z x f θ yz f x y f θ x z f θ y z f θ θ 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π Chapter 3 – Numerical tools for soil-structure interaction 104 Figure 3.26: Tunnel compliances for the layered half-space. Solid lines = 2.5D BEM (real part – blue; imaginary part – red). Black circles = FEM 0 1 2 3 -0.1 -0.05 0 0.05 0.1 0 1 2 3 -0.1 -0.05 0 0.05 0.1 0 1 2 3 -0.1 -0.05 0 0.05 0.1 0 1 2 3 -0.4 -0.2 0 0.2 0.4 0 1 2 3 -0.2 -0.1 0 0.1 0.2 0 1 2 3 -0.4 -0.2 0 0.2 0.4 0 1 2 3 -0.05 0 0.05 0 1 2 3 -0.1 -0.05 0 0.05 0.1 0 1 2 3 -0.02 -0.01 0 0.01 0.02 0 1 2 3 -0.05 0 0.05 0.1 0.15 0 1 2 3 -0.04 -0.02 0 0.02 0.04 0.06 0 1 2 3 -0.06 -0.04 -0.02 0 0.02 0.04 xx f 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π 2 y k π zz f z z f θ θ y x f θ z x f θ yz f x z f θ y z f θ θ yy f y y f θ θ x y f θ x x f θ θ Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 105 3.5.2 Example 2 – slab free in space For the second example an infinitely long slab with width 1m L = and thickness 0.1m H = that is free in space is considered. The material properties of the slab are: mass density 3 1kg m ρ = ; shear modulus 1Pa G = ; Poison’s ratio 0.25 ν = ; material damping 0.01 ξ = . The slab is submitted to a vertical line load at its top (at the middle alignment) and the displacements of the top right corner, center right point and bottom right corner are calculated. The excitation frequency is 2 rad/s ω π = and the range of wavenumbers of the line load is [0,6 ] y k π = rad/m . The displacements of the points referred to above are calculated using a 2.5D FEM procedure and using a coupled 2.5D BEM-FEM procedure, and then compared. For the first approach the cross section of the slab is divided into a regular mesh of 20 8 × solid elements of quadratic expansion (8 nodes per element), while for the second approach the cross section of the slab is divided into two sub-domains: a FEM sub-domain, with dimensions 1m 0.05m × and divided into a regular mesh of 20 4 × solid elements of quadratic expansion; and a BEM sub-domain with the same dimensions, whose lateral boundaries are divided into 4 boundary elements each and whose interface between the BEM and the FEM sub-domains is divided into 20 boundary elements. The boundary elements are of quadratic expansion (3 nodes per element). The fundamental solutions of an elastic layer free in space, calculated with the TLM (the elastic layer is divided into 32 thin-layers of quadratic expansion), are used to nurture the boundary elements: this choice for the fundamental solutions avoids the discretization of the free lower surface of the slab, thus reducing the cost of computation of the BEM matrices. The results obtained with the two approaches are compared in Figure 3.27-3.29 for the 3 points considered. As can be observed, the agreement is very good even for the center right point, which belongs to the BEM-FEM interface, and for the bottom left corner, which belongs to the BEM domain. The good quality of the results validates the two procedures. As a concluding remark regarding this example, it is important to be aware that the use of the 2.5D BEM-FEM is inefficient when compared to the 2.5D FEM. The structure under analysis is finite and of relatively small dimensions, which means that the reduction in the number of degrees of freedom achieved with the BEM does not compensate the computational cost associated with the calculation of the BEM matrices. This example is considered herein solely for validation purposes, and not to demonstrate the advantages of the 2.5D BEM. Chapter 3 – Numerical tools for soil-structure interaction 106 Figure 3.27: Displacements of the top right corner of the slab: a) u x ; b) u y ; c) u z . Solid lines = 2.5D BEM (real part – blue; imaginary part – red). Black circles = FEM Figure 3.28: Displacements of the center right point of the slab: a) u x ; b) u y ; c) u z . Solid lines = 2.5D BEM (real part – blue; imaginary part – red). Black circles = FEM 0 5 10 15 20 -0.2 -0.1 0 0.1 0.2 0.3 k y u x a) 0 5 10 15 20 -1.5 -1 -0.5 0 0.5 1 k y u y b) 0 5 10 15 20 -10 -5 0 5 10 15 k y u z c) 0 5 10 15 20 -6 -4 -2 0 2 4 k y u x a) 0 5 10 15 20 -4 -2 0 2 4 k y u y b) 0 5 10 15 20 -10 -5 0 5 10 15 k y u z c) Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 107 Figure 3.29: Displacements of the bottom right corner of the slab: a) u x ; b) u y ; c) u z . Solid lines = 2.5D BEM (real part – blue; imaginary part – red). Black circles = FEM 3.6 Conclusions In this chapter, the numerical tools used for the solution of soil-structure interaction problems are presented and validated. To solve the interaction between the track and the soil, a coupled 2.5D BEM-FEM procedure is used, and for the case of the interaction between the building and the soil, a coupled 3D BEM-FEM procedure is used. Both procedures are based on the fundamental solutions obtained with the TLM, which is described in chapter 2. The 2.5D BEM-FEM procedures are developed in the space-wavenumber-frequency ( , y k ω ) domain. In order to use the responses obtained with this methodology in a coupled 3D BEMFEM procedure, the responses must first be transformed to the space-frequency ( , y ω ) domain, which is accomplished by an inverse Fourier transform. After having calculated the responses in the ( , y k ω ) domain for a discrete sample of the wavenumber y k , the inverse transform can be calculated numerically by means of a summation. It is important to recall that when the proposed 2.5D BEM methodology is used, the coefficients for the BEM matrices can be calculated in closed form expressions. This fact leads to fast and precise calculations of such coefficients. The drawback of this approach is the time required to calculate the eigenmodes of the soil, which can become large when the fundamental solutions are needed at deep positions. Nevertheless, for each soil profile, the eigenmodes only have to be calculated once for each frequency. Then, they can be stored and reused to analyze different configurations of tracks, buildings and countermeasures. Concerning the BEM procedures, the step that consumes more time is the calculation of the BEM matrices P and U . The components of these matrices are obtained by applying a load 0 5 10 15 20 -4 -2 0 2 4 6 k y u x 0 5 10 15 20 -4 -2 0 2 4 k y u y 0 5 10 15 20 -10 -5 0 5 10 15 k y u z Chapter 4 – Invariant structures subjected to moving loads and moving vehicles 114 When the load speed is near cr V and the damping is neglected, the roots j ω are very close to the real axis and have approximately the same value (apart from the sign) and so the response propagates both ahead and behind the load with roughly the same wavelength, which is approximately cr re 2V π ω , being ( ) 2 re 2 V m EI ω ≈ . The presence of damping makes the wavelengths ahead and behind the load different, but their values are still close to each other, as can be inferred from Figure 4.2b-c. Finally, for cr V V > and again for null damping, there are two distinct pairs of real roots: 1 ω ± and 2 ω ± ( 1 2 ω ω > ), being 1 ω associated with the waves ahead of the load and 2 ω associated with the waves behind the load position, thus turning the wavelength ahead of the load smaller than the wavelength behind the load. As V increases, 1 ω increases and 2 ω decreases (tending to a minimum of k m ) and consequently the wavelength ahead of the load shortens, while the wavelength behind the load increases. In the case under study, the presence of damping makes the roots j ω complex, but the features explained in the last sentence are still present, as can be observed in Figure 4.2d. Also observable in Figure 4.2a-d is the amplification of the displacements when the load speed approaches cr V . That fact is confirmed in Figure 4.3, where the maximum beam displacement is plotted as a function of the load speed V for two distinct damping scenarios: 3 30 10 Ns/m c= × (blue line) and 0 Ns/m c = (red line). The effect of the damping is similar to the effect of a dashpot on one-degree-of-freedom systems, which reduces the amplification of the response at the resonance frequency of the systems. Figure 4.3: Beam maximum displacements as a function of the load speed. Blue line = with damping; red line = no damping 4.2.3 Oscillating moving loads Consider now that a load moving with constant speed V and with time varying amplitude defined by ( ) ( ) 0 0 0 , exp i f t A t ω ω = crosses the longitudinal section 0 y y = at the instant 0 t = 0 200 400 600 800 1000 0 0.002 0.004 0.006 0.008 0.01 V [m/s] max( u) Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 115 (the forcing frequency 0 ω is added as an argument for convenience). The function ( ) 0 , , p y t ω associated with such load is ( ) ( ) 0 i 0 0 0 , , e t p y t A y Vt y ω ω δ = − − (4.15) and its wavenumber-frequency content ( ) 0 , , y p k ω ω ɶ is ( ) ( ) ( ) ( ) 0 0 0 0 i i-i 0 0 0 -i i i 0 0 0 , , e e d e d e e d 2 e y y y y k y tt y t k V k y k y y p k A y Vt y y t A t A k V ωω ω ω ω ω δ π δ ω ω +∞ +∞ −∞ −∞ +∞ − − −∞ = − − = = = + − ∫ ∫ ∫ ɶ (4.16) The insertion of (4.16) and (4.3) in equation (4.1) yields ( ) ( ) ( ) ( ) 0 ii 0 0 0 1 , , , e d e d 2 y k y y t y y y h y t A h k k V k ω ω ω δ ω ω ω π +∞ +∞ − − −∞ −∞ = + − ∫ ∫ ɶ (4.17) and after solving equation (4.17) for the inner integral, the following expression is obtained ( ) ( ) 0 0 0 0 0 0 0 , , i -i i 0 0 0 i i 0 0 1 , , e , e e d 2 e , e d 2 h y y y y y t V V y y y y tV V A h y t h V V Ah V V ω ω ω ω ω ω ω ω ω ω ω ω π ω ω ω ω π − − +∞ −∞ −   −+∞ −     −∞ −   = =     −   =     ∫ ∫  ɶ ɶ (4.18) The force ( ) 0 , f t ω described above has no physical meaning since it contains an imaginary component, which results from the exponential factor ( ) 0 exp i t ω . However, this type of expression can be used to define real functions of the type cosine or sine, as stated in the next equations ( ) ( ) ( ) ( ) 0 0 0 0 i -i C 0 0 0 C i -i S 0 0 0 S e e , cos 2 e e , sin 2i t t t t f t A t A f t A t A ω ω ω ω ω ω ω ω + = = − = = (4.19) Hence, in order to obtain the response ( ) C 0 , , h y t ω for cosine type loads or the response ( ) S 0 , , h y t ω for sine type loads, equation (4.18) must first be used to calculate ( ) 0 , , h y t ω and afterwards, since ( ) ( ) 0 0 , , , , h y t h y t ω ω − = , the real or imaginary component must be retained according to the following ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 0 0 C 0 0 0 0 S 0 0 , , , , , , Re , , 2 , , , , , , Im , , 2i h y t h y t h y t h y t h y t h y t h y t h y t ω ω ω ω ω ω ω ω + − = =     − − = =     (4.20) Similarly to the case of constant moving loads, for oscillating moving loads no integral needs to be evaluated in order to obtain the space-frequency domain response. Instead, the response in that domain is obtained simply by evaluating the transfer function ( ) , y h k ω ɶ for the Chapter 4 – Invariant structures subjected to moving loads and moving vehicles 116 wavenumber-frequency pair ( 0 ( ) , V ω ω ω −). In addition, the response of points that move with the same speed as the load is given by ( ) PP 0 ii 0 0 0 P , e , e d 2 yy tVV A h y Vt y y t h V V ωω ω ω ω ω π   +∞ −     −∞ −   = + − =     ∫ ɶ (4.21) in which P 0 y Vt y y = + − is the distance between the moving point and the source. The previous equation clearly shows that the response fields move together with the load and oscillate with frequency 0 ω . Beam on a Kelvin foundation subjected to a cosine type moving load For a load ( ) ( ) 0 0 cos f t A t ω =, the beam displacements are calculated with ( ) ( ) C 0 0 , , Re , ,u y t u y t ω ω =    (4.22) in which ( ) ( ) 0 0 0 i 3i 0 044 4 2 4 0 e , , e d 2i y y t y y V V AV u y t EI kV cV mV ω ω ω ω πω ω ω ω −   −   −+∞   −∞ =− + + − ∫ (4.23) Since the integrand in (4.23) is a regular expression and since it is bounded in the complex plane, the integral can be evaluated by means of contour integration, resulting in ( ) ( ) ( ) 0 0 0 0 34 4i i 0 0 011 imag 0 1 , , i e sign e j j y y y y tV V jkj k k j y y tV AV y y u y t t EI V ω ω ω ωω ω −   −−     == ≠ −   − >       −     = −     −       ∑ ∏ (4.24) In the previous equation, the values j ω represent the roots of the polynomial ( ) 44 4 2 4 0 i 0 j j j EI kV cV mV ω ω ω ω − + + − = (4.25) Similarly to the case of constant moving loads, when damping is null, there is a velocity above which the polynomial (4.25) presents real roots. However, in this case, for low frequencies 0 ω , the minimum velocity that yields real roots only yields two of the poles real, and so there is a second velocity above which all the roots are real. At these two critical velocities, two roots of (4.25) become real and with the same value (double roots), which results in the amplification of the displacements of the beam. In Figure 4.4, the (logarithm of) maximum of the displacements C u is plotted as a function of the excitation frequency 0 ω and load speed V for a beam on a Kelvin foundation with the same properties of the system described in section 4.2.2 (with 0 Ns/m c = ). Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 117 Figure 4.4: Logarithm of the beam maximum displacement as a function of the load speed V and excitation frequency 0 ω ( 0 Ns/m c = ) For 0 0 ω = , there is only one critical velocity which is calculated with (4.14). For, 0 0 ω ≠ , to determine the load speeds is cumbersome, but by rearranging equation (4.25) it can be concluded that these velocities correspond to the local real minima of ( ) 4 0 4 2 , EI k V m k m ω ω ω ω − = > − (4.26) When 0 k m ω >, there is only one critical velocity. Next, the beam displacements are calculated considering the excitation frequency 0 100 rad/s ω = and considering the load speeds 300 V = , 500 V = and 700m/s V = . The properties of the beam and foundation are the same of the example of section 4.2.2, and it is assumed that the load crosses the cross-section 0 y = at the instant 0 t = ( 0 0 y = ). The results obtained with the time domain FEM procedure are not shown, but is was observed that the agreement between the FEM results and equations (4.22) and (4.23) is very good. The displacements at 0 t = are plotted in Figure 4.5. For load speeds below the lower critical velocity (Figure 4.5a), the displacements evanesce with the distance to the source and the decay rate of the displacements depends on the proximity of the load speed to the critical velocity. For load speeds between the two critical velocities (Figure 4.5b), the displacements propagate both ahead ( 0 y > ) and behind ( 0 y < ) the load, being the wavelength of the displacements shorter ahead than behind. For load speeds greater than the highest critical velocity (Figure 4.5c), the displacements propagate mostly behind the load and the wavelengths are shorter ahead than behind. Note that there are two major wavelengths for each side of the response, a consequence of the existence of four distinct real roots in the polynomial (4.25) (for non-oscillating loads, the roots are paired up in groups of 2, 1 ω ± and 2 ω ± , and so the two distinct wavelengths are not observed). Chapter 4 – Invariant structures subjected to moving loads and moving vehicles 118 Figure 4.5: Beam displacements induced by harmonic moving loads In Figure 4.6, the maximum displacements observed at the beam are plotted as a function of the load speed for the conditions referred to above (blue line) and for null damping (red line). The two critical velocities can be identified for the case of no damping, while for the case of damping, the critical velocities appear to merge into one. Such feature is observed for the damped case because the critical velocities are near each other and because the large amount of damping used not only reduces the maximum displacements as it also widens the “bell shaped response” associated with each velocity. If the damping is reduced or if the excitation frequency 0 ω is increased (and consequently the critical velocities are moved farther apart), then, even for the damped case, the two critical velocities can be noticed. The example solved in this subsection considers cosine type loads. For sine type loads, the same conclusions can be drawn and the response of the structure is practically the same, being it simply shifted in time and in space. In the next subsection, structures for which the response fields cannot be determined analytically are studied. -20 -10 0 10 20 -5 0 5 10 15 20x 10 -4 y u C (y,t) a) V = 300 [m/s] -20 -10 0 10 20 -2 -1 0 1 2 3 x 10 -3 y u C (y,t) b) V = 500 [m/s] -30 -20 -10 0 10 20 30 -2 -1 0 1 2 x 10 -3 y u C (y,t) c) V = 700 [m/s] Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 119 Figure 4.6: Beam maximum displacements as function of the load speed. Blue line – with damping; Red line – no damping 4.2.4 Examples Next, equations (4.9) and (4.18)-(4.20) are used to calculate the response of longitudinally invariant structures of which the wavenumber transfer functions (and therefore the integrands of the mentioned equations) cannot be determined analytically. The case study considered in this subsection consists in a flexible slab resting on the surface of a half-space (homogeneous or layered) that is subjected to a vertical load ( ) f t that moves with speed V along the middle top alignment (point A), as exemplified in Figure 4.7 (it is assumed that the load crosses the section 0 y = at 0 t = , i.e., 0 0 y = ). The material properties of the slab are: density 3 2145 kg/m ρ =; Young’s modulus 30 GPa E = ; and Poisson’s ratio 0.2 ν = . The slab is modeled with 8 four-node volume elements with dimensions 2 0.25 0.3 m × and the interface between the slab and the half-space is divided into 8 boundary elements of constant expansion. Figure 4.7: Slab resting on a half-space submitted to a moving load Example 1 – Slab on a homogeneous half-space subjected to a constant moving load In this first example, the slab is subjected to a constant moving load defined by ( ) 1000 N f t = and the foundation, which consists in a homogeneous half-space, is given the following properties: density 3 1800 kg/m ρ =; shear modulus 0.1125 GPa G = ; Poisson’s ratio 0.25 ν = 1m 1m 0.3m A B ( ) f t V C 1m D E 4m 5m x y z 0 200 400 600 800 1000 0 0.002 0.004 0.006 0.008 0.01 V [m/s] max( u C ) Chapter 4 – Invariant structures subjected to moving loads and moving vehicles 120 (the corresponding body wave velocities are s 250 m/s C= and p 433 m/s C=). A small amount of hysteretic damping p s 0.02 ξ ξ = = is considered, which renders complex wave velocities ( ) p p p 1 sign 2i C C ω ξ = + and ( ) s s s 1 sign 2i C C ω ξ = + (Dominguez, 1993). For constant moving loads, the time domain response fields ( ) , h y t are obtained through the evaluation of eq. (4.9), but since in this example the transfer functions ( ) , y i h k ω ɶ cannot be determined analytically, then the integral in (4.9) has to be evaluated numerically, being approximated with ( ) ( ) 0 i 0 , , e 2 i y y NtV i i i N A h y t h V V ω ω ω ωω ω π −   −     =− ∆ ≈ ∑ ɶ (4.27) Due to the conjugate property ( ) ( ) , , y i y i h k h k ω ω − − = ɶ ɶ , eq. (4.27) can be further simplified to ( ) ( ) ( ) ( ) ( ) 00 0 1 , Re , cos Im , sin Ny y y y t t i i i i i i V V i A h y t h V h V V ω ω ω ωω ω ω ω π − −     − −             = ∆       ≈ −         ∑ ɶ ɶ (4.28) In the following, equation (4.28) is used to calculate the time domain displacements at the cross-section 0 y = and for 100 m/s V = at point B, situated at the edge of the slab, and at points C, D and E, situated at the surface of the half-space. Two TLM models are tested: in the first model (TLM-1), the half-space is modeled with just one PML with parameters 2 η = , 4 Ω = , 2 m = and 10 N = (see chapter 2); in the second model (TLM-2), the half-space is modeled with an elastic layer of thickness 2 m H = , divided into 40 quadratic thin-layers, and with a PML (same parameters as for model TLM-1). The displacements obtained with these TLM models are compared with the displacements obtained using a time domain methodology (TD) (dos Santos et al., 2010a; dos Santos et al., 2010b) and with the displacements obtained using eq. (4.28) together with the transfer functions ( ) , y i h k ω ɶ obtained from a 2.5D BEM-FEM procedure based on the stiffness matrices of Kausel and Roesset (1981) (SM). For the procedures based on eq. (4.28), 1500 frequencies with a step of 0.1 Hz are used. For the time domain procedure, a 100m long 3D model divided into 400 longitudinal sections and a time step of 0.002 s are used. Figure 4.8 plots the transverse displacements of the points B, C, D and E obtained with the four mentioned approaches, Figure 4.9 plots the longitudinal displacements and Figure 4.10 plots the vertical displacements. From Figure 4.8, it can be observed that the four approaches yield significantly different transverse displacements: at the slab (Point B), the TLM-1, the TLM-2 and the SM results tend approximately to the same maximum value, but the shapes of the curves are different; still regarding point B, the TD solution differs both in shape and sign from the remaining solutions (a justification for this could not be found); at the surface of the half-space (points C, D and E), the results of the SM model tend the follow the TD results while the TLM-1 and the TLM-2 yield smaller displacements. Despite this fact, the passage of the load at the crosssection is noticed in all four approaches. As for the longitudinal and vertical displacements represented in Figure 4.9 and Figure 4.10, respectively, a better agreement is, in general, observed. Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 121 Figure 4.8: Transverse (x) displacements: blue = TLM-1; red = TLM-2; black = TD; green = SM Figure 4.9: Longitudinal (y) displacements: blue = TLM-1; red = TLM-2; black = TD; green = SM -0.2 -0.1 0 0.1 0.2 0.3 -1 -0.5 0 0.5 1x 10 -7 t [s] u y [m] - 0.2 -0.1 0 0.1 0.2 0.3 -5 0 5x 10 -8 t [s] u y [m] -0.4 - 0.2 0 0.2 0.4 -4 -2 0 2 4x 10 -8 t [s] u y [m] -0.5 0 0.5 -4 -2 0 2 4x 10 -8 t [s] u y [m] Point B Point C Point D Point E - 0.2 - 0.1 0 0.1 0.2 0.3 -6 -4 -2 0 2x 10 -8 t [s] u x [m] -0.2 -0.1 0 0.1 0.2 0.3 -15 -10 -5 0 5x 10 -8 t [s] u x [m] -0.4 -0.2 0 0.2 0.4 -6 -4 -2 0 2x 10 -8 t [s] u x [m] -0.5 0 0.5 -4 -3 -2 -1 0 1x 10 -8 t [s] u x [m] Point B Point C Point D Point E Chapter 4 – Invariant structures subjected to moving loads and moving vehicles 122 Figure 4.10: Vertical (z) displacements: blue = TLM-1; red = TLM-2; black = TD; green = SM When comparing the TD results (black lines) with the remaining approaches, it is noticed that for the earlier moments ( 0.2 t < s) the responses present very distinct behaviors. These differences are justified by the finite length of the 3D model used in the TD approach, a characteristic that violates the assumption of invariant cross-section considered in the 2.5D models. In this way, while in the TD approach the entrance of the load in the finite element model induces a transient phenomenon, in the 2.5D approaches, since it is assumed that the load travels from minus infinity to plus infinity, the phenomenon is not present. The transient phenomenon dissipates due to internal damping and therefore after some time its contribution is minimal. Apart from the transient phenomenon, the remaining differences may be justified by the longitudinal and temporal discretizations required by the TD approach: the TD and the SM approaches are based on the fundamental solutions obtained with the stiffness matrices and therefore, theoretically, they should yield the same values. This hypothesis has not been tested because to make the 3D mesh longer in the longitudinal direction and/or with thinner elements and smaller time-steps renders the calculation unfeasible (for the current 100 m long model, more than two days were required to calculate the passage of the load from one edge of the model to the other edge). The differences between the results of the TLM-1, TLM-2, and SM approach can only be justified by differences in the transfer functions ( ) , y h k ω ɶ . Figure 4.11 compares the transfer functions ( ) / , x u V ω ω ɶ calculated with the three approaches at the four points and, as expected, they present discrepancies, which justify the distinct curves in Figure 4.8. Furthermore, when comparing the blue and red lines (TLM-1 and TLM-2) with the green line (SM), one realizes that the differences are concentrated in the lower frequency range ( 20 rad/s 3.2 Hz f ω < → < ), -0.2 -0.1 0 0.1 0.2 0.3 -15 -10 -5 0 5x 10 -7 t [s] u z [m] -0.2 - 0.1 0 0.1 0.2 0.3 -6 -4 -2 0 2x 10 - 7 t [s] u z [m] -0.4 -0.2 0 0.2 0.4 -20 -15 -10 -5 0 5x 10 -8 t [s] u z [m] - 0.5 0 0.5 -15 -10 -5 0 5x 10 - 8 t [s] u z [m] Point B Point C Point D Point E Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 123 which leads to the conclusion that these differences are caused by the inappropriate simulation of the infinite domain at the low frequency range. The comparison of the longitudinal and vertical responses is not shown here, but it was observed that the differences are smaller than for the transverse component, which justifies the better agreement obtained. Figure 4.11: Transverse (x) displacements in the wavenumber-frequency domain: solid line = real part; dashed line = imaginary part; blue = TLM-1; red = TLM-2; green = SM Note that for the low frequencies, the PML used to model the half-space is divided into thinlayers that are thicker than the widths of the slab and of the boundary elements (the PML and the thin-layers are made proportional to the shear wavelength) and so the boundary integrals calculated as explained in chapter 3 do not yield accurate results. Also, at low frequencies, when the waves, that are mostly evanescent, enter the PML, they propagate instead of evanesce, and that causes the waves to reflect at the bottom of the PML and return to the surface, which contaminates the results. This is also the reason why the TLM-2 model yields better results than the TLM-1 model: the TLM-2 model contains an elastic layer divided into thinner thin-layers where the waves evanesce before they reach the PML. A third TLM model in which the elastic part is made thicker (10 m divided into 100 quadratic thin-layers) has been tested and in this case the transfer functions ( ) / , x u V ω ω ɶ (and remaining components) approximate more closely the green line. Hence, the problem of simulating the half-spaces for low frequencies can be solved by using a sufficiently thick elastic thin-layer on top of the PML. An alternative approach to overcome the above mentioned differences would be to replace the PML by the static half-space matrices derived by Kausel and Seale (1987), thus neglecting the mass of the half-space (in the very low frequency range, this simplification should not introduce significant errors). However, the terms derived in the mentioned work change the 0 100 200 300 400 -15 -10 -5 0 5x 10 -8 ω [rad/s] u x (ω/V, ω) [m] 0 100 200 300 400 -8 -6 -4 -2 0 2x 10 -7 0 100 200 300 400 -10 -5 0 5x 10 - 7 0 100 200 300 400 -10 -5 0 5x 10 -7 u x (ω/V, ω) [m] u x (ω/V, ω) [m] u x (ω/V, ω) [m] ω [rad/s] ω [rad/s] ω [rad/s] Point B Point C Point D Point E Chapter 4 – Invariant structures subjected to moving loads and moving vehicles 130 Figure 4.19: Maximum vertical displacement of point B as a function of the load speed V for: 0 0 ω = (gray); 0 10 rad/s ω π = (blue); 0 20 rad/s ω π = (red); 0 30 rad/s ω π = (black); and 0 40 rad/s ω π = (green) Figure 4.20: Vertical displacements for 100 m/s V = ( 0 40 rad/s ω π =): the displacements are multiplied by the shear modulus of the half-space (values in N/m) 100 150 200 250 300 350 400 450 500 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 x 10 - 6 V [m/s] max( u z ) [m] Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 131 Figure 4.21: Vertical displacements for 300 m/s V = ( 0 40 rad/s ω π =): the displacements are multiplied by the shear modulus of the half-space (values in N/m) Figure 4.22: Vertical displacements for 500 m/s V = ( 0 40 rad/s ω π =): the displacements are multiplied by the shear modulus of the half-space (values in N/m) Chapter 4 – Invariant structures subjected to moving loads and moving vehicles 132 For 100 m/s V = (Figure 4.20) it can be observed that the waves propagate away from the load and that the wave-fronts present an elliptical shape, being the wavelengths ahead of the load shorter than the wavelengths behind the load. Such aspect is characteristic of waves induced by loads moving with speeds below the surface waves of the medium. Unlike the case of a constant moving load with the same speed (Figure 4.13), in this example the values of the displacements alternate between positive and negative, a consequence of the oscillatory nature of the load. Also, the reduction of displacements with the distance to the source is smaller, which emphasizes the importance of the dynamic loads in the calculation of wave-fields at remote positions. For 300 m/s V = and for 500 m/s V = (Figure 4.21 and Figure 4.22, respectively) the wavefronts form a cone, which is also observed for the case of constant moving loads, and that is a consequence of the load moving faster than the surface waves of the foundation. Similarly to constant loads, the cone is shaper for 500 m/s V = than for 300 m/s V = , and for 500 m/s V = there is a second cone that results from the fact that the load moves faster than the pressure wave of the foundation (this second cone cannot be noticed in Figure 4.22). Example 3 – Critical speeds for layered foundations The first two examples considered homogeneous foundations. In this example, the critical load speeds are calculated for different configurations of the foundation, namely: a) a stiffer layer on top of a softer half-space ( lay half 0.225 GPa, 0.1125 GPa G G = = ); b) a softer layer on top of a stiffer half-space ( lay half 0.05625 GPa, 0.1125 GPa G G = = ). The thickness of the upper layers is 2 m H = (divided into 40 quadratic thin-layers), and the material properties are the same of the homogeneous half-space used in examples 1 and 2 (except for the shear modulus, which is as indicated in the previous sentence). According to Dieterman and Metrikine (1997), resonance occurs when the load speed equals the group velocity of the waves generated by the load, which in other words means that resonance takes place when the integrating path 0 y Vk ω ω = + or 0 y Vk ω ω = − + is tangent to a dispersion curve of the foundation. The dispersion curves of a layered domain represent its free vibration modes and are defined by the curves ( , ( ) k k ω ) that yield the system singular, i.e., undetermined. These curves are associated with the steady waves that are originated due to the existence of boundary conditions (such as free surface or interface between layers). These waves propagate horizontally with phase velocity ph ( ) V k k ω = and group velocity gr ( ) V k k ω = ∂ ∂ . For the case of a homogeneous full-space, there are no such steady waves because the only existent boundary condition is the radiation of waves to infinity. As for the case of a homogeneous half-space, the existence of a free surface originates the well-known Rayleigh wave, whose phase and group velocities are constant (Rayleigh wave speed). Under this scenario, the system is classified as non-dispersive. For the cases of layered systems, the number of existing waves and their phase velocity depend on the frequency being considered. The system is then said to be dispersive. The Love waves, which are observed in anti-plane systems consisting of a layer on a half-space, and the Stonely waves, which can exist at the interface between two adjacent layers for a certain range of frequencies, are examples of such dispersive waves (Erigen and Suhubi, 1975). Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 133 Besides being useful for the calculation of critical load speeds, the dispersion curves of a foundation can also be used to investigate some properties of the response of the layered domain, especially for remote positions where the evanescent waves are not noticed. For instance, for the case of non-moving harmonic loads with excitation frequency 0 ω , the waviness of the response at remote positions is characterized by the wavelengths 2 j k π , being j k the wavenumbers at which the line 0 ω ω = intercepts the dispersion curves. For moving loads, the dispersion curves indicate which frequencies dominate the response at the free-field, which correspond to the frequencies at which the lines 0 y Vk ω ω = + and 0 y Vk ω ω = − + intercept the dispersion curves. Using the TLM, the dispersion curves are obtained through the calculation of the real modes k of the SH and SVP eigenvalue problems as a function of the frequency ω . However, when in the presence of PMLs to simulate half-spaces, the dispersion curves obtained with this method may not have any physical meaning, since the PML is not more than an artifact to simulate a half-space up to a certain distance. Alternatively, making use of the stiffness matrices of Kausel and Roesset (1981), the pairs ( , ( ) k k ω ) are found by forcing the assembled stiffness matrices to become singular, or in other words, by forcing the determinant of the matrix to be null. This last procedure has been used in this work together with search techniques to calculate the dispersion curves of the layered system considered as foundations. The obtained curves are represented in Figure 4.23 together with the corresponding phase velocities ph V and group velocities gr V . Figure 4.23 also plots the expected critical load speeds that, according to Dieterman and Metrikine (1997), are calculated by ( ) ( ) 0 gr k kV k ω ω = − (4.31) The first feature that can be observed in Figure 4.23 is that for layered half-spaces no wave propagates with phase velocity higher than the shear wave velocity of the half-space ( s 250 m/s C =). That is so because for such waves to keep propagating, energy has to come from inside the half-space (otherwise the waves would propagate faster than it is admissible in the half-space), and that violates the radiation condition. Due to the reason explained above, the inversely dispersive domain a) only presents one real pole, and that pole becomes complex for frequencies 52 ω π > (rad/s), which is approximately the frequency at which the phase velocity of the surface wave reaches the shear wave velocity of the half-space. In addition, the lowest critical load speed of this domain occurs for 0 0 ω = and equals the Rayleigh wave velocity of the half-space, as can be inferred from the last row of Figure 4.23. The highest admissible critical load speed is approximately cr 263 m/s V ≈, which occurs for a excitation frequency 0 2.7 ω π ≈ rad/s. Above that excitation frequency, there are no critical load speeds. As for the normally dispersive domain b), since the layer is softer than the half-space, there is a vast set of waves that can coexist at the interface between the layer and the half-space. Since one of the lower poles stabilizes at the phase velocity ph 162 m/s V ≈ (Rayleigh wave speed of the layer), for the frequency of excitation 0 0 ω = , such velocity is one of the critical load speeds of the system. Nevertheless, that is not the minimal admissible critical load speed, as can be inferred from the last row of Figure 4.23: the lowest mode admits a critical speed of cr 142 m/s V ≈ that occurs for the excitation frequency 0 17 ω π ≈. In practice, this must be the Chapter 4 – Invariant structures subjected to moving loads and moving vehicles 134 maximum speed for the vehicle to circulate without causing the resonance of the track-soil system. Figure 4.23: Dispersion curves, phase velocities, group velocities and expected critical load speeds for the layered domain a) (left) and for the layered domain b) (right) 0 0.2 0.4 0.6 0.8 0 20 40 60 Dispersion Curves k y ω/π 0 20 40 60 230 235 240 245 250 Phase Velocity ω/π V ph [m/s] 024 6 8 0 100 200 300 400 k y ω/π 0 100 200 300 400 160 180 200 220 240 260 ω/π V ph [m/s] 0 20 40 60 230 240 250 260 270 Group Velocity ω/π V gr [m/s] 012 3 230 240 250 260 270 Critical Load Speed ω 0 /π V cr [m/s] 0 100 200 300 400 100 150 200 250 ω/π V gr [m/s] 0 50 100 150 100 150 200 250 ω 0 /π V cr [m/s] Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 135 Next, the critical load speeds are calculated using equations (4.28) and (4.30). With that intention, the maximum displacements of point B are calculated for a set of values of 0 ω and for the load speeds ranging from 80 to 480 m/s. The considered values of 0 ω are: 0 0 ω = , 0 rad/s ω π =, 0 2 rad/s ω π = and 0 4 rad/s ω π = for the layered domain a); 0 0 ω = , 0 10 rad/s ω π =, 0 15 [rad/s] ω π = and 0 40 rad/s ω π = for the layered domain b). The results are shown in Figure 4.24. Figure 4.24: Maximum vertical displacement of point B as a function of the load speed V and excitation frequency 0 ω for the layered domain a) and layered domain b) The analysis of the displacement envelope of the layered domain a) (left image of Figure 4.24) indicates that for 0 0 ω = , 0 rad/s ω π = and 0 2 rad/s ω π = the critical load speeds are effectively those predicted in Figure 4.23. However, there is no evidence of the existence of the higher critical speed associated with the upper part of the blue line of Figure 4.23 (bottom left figure) for 0 rad/s ω π = and for 0 2 rad/s ω π =. Also, there is a peak at 270 m/s V ≈ for 0 4 rad/s ω π = that according to Figure 4.23 does not correspond to a critical load speed. As for the layered domain b), for the excitation frequencies 0 10 rad/s ω π = and 0 15 rad/s ω π =, the critical load speeds cr 160 m/s V ≈ and cr 150 m/s V ≈ associated with the upper part of the blue curve of Figure 4.23 (bottom right figure) can be identified in the right image of Figure 4.24. Nevertheless, the critical load speeds associated with the lower part of that curve cannot be distinguished. Also, for these excitation frequencies, the amplification is greatest for the critical load speeds associated with the red line (around 240 m/s). For 0 0 ω = , it can be observed that the critical load speed associated with the upper part of the blue curve is slightly shifted to the left, being its value cr 215 m/s V ≈ instead of the expected 230 m/s . A local maximum can also be observed at the load speed 250 m/s V = , which corresponds to the group velocity with which the second and higher poles become real, and equals the shear wavevelocity of the half-space. Finally, for the excitation frequency 0 40 rad/s ω π = no critical load speed can be detected in Figure 4.24. The differences between the critical load speeds calculated according to Dieterman and Metrikine (1997) and the critical speeds detected in Figure 4.24 can be justified by the presence of the slab, since the slab, which is not considered in the mentioned work, changes the behavior of the foundation, namely for the higher frequencies (Steenbergen and Metrikine, 0 100 200 300 400 500 0.4 0.6 0.8 1 1.2 1.4 1.6 x 10 - 6 V [m/s] max( u z ) [m] Layered domain a) ω 0 = 0 ω 0 = π ω 0 = 2 π ω 0 = 4 π 0 100 200 300 400 500 1 1.5 2 2.5 3 3.5 4 x 10 -6 V [m/s] max( u z ) [m] Layered domain b) ω 0 = 0 ω 0 = 10 π ω 0 = 15 π ω 0 = 40 π Chapter 4 – Invariant structures subjected to moving loads and moving vehicles 136 2007) (note that when the integrating paths 0 cr V k ω ω = + and 0 cr V k ω ω = − + are tangent to the dispersion curves at low frequencies ω , the displacement envelope presents a local maximum for cr V V = ). Nonetheless, despite the differences, the critical load speeds calculated according to Dieterman and Metrikine (1997) can be very helpful to estimate the vehicle speeds that can cause the resonance of the track-soil system. 4.2.5 Conclusions In this section, the equations used to calculate the response of invariant structures subjected to loads moving with constant speed are derived and discretized so that they can be used in cases in which the transfer functions cannot be determined in closed-form expressions. Important to retain is that in order to transform the transfer functions from the wavenumber-frequency domain to the space-frequency domain, no integral needs to be solved. Only if one attempts to obtain the space-time domain response it is needed to solve an integral. Theoretical examples consisting of a beam on a Kelvin foundation are used to validate the equations and some conclusions concerning the critical load speeds are presented. Then, the example of a slab resting on the surface of a half-space is considered and the results obtained with the proposed methodology are compared with the results obtained with a time domain 3D FEM procedure. It is observed that the proposed procedure yields slightly different results for moving loads with constant magnitude (mostly due to the incorrect simulation of the halfspace at very low frequencies) and that it yields excellent results for oscillating moving loads. This dynamic component of the load is crucial for a correct prediction of the vibrations at remote positions. The influence of the speed of the load on the response of the slab-foundation system is also investigated, and it is concluded that while for foundations consisting of homogeneous halfspaces there is a critical load speed that coincides with the Rayleigh wave velocity of the halfspace (but only if the load magnitude is constant in time), for foundations consisting of layered half-spaces the critical load speeds depend on the stratification of the domain (in the sense that it defines the dispersion curves) and on the frequency of the moving load. Hence, both the stratification of the foundation and the dynamic component of the loads are important aspects to consider in the calculation of the response of structures subjected to moving loads or vehicles. 4.3 Moving vehicles 4.3.1 Introduction The forces that a moving vehicle transmits to the supporting structure can be divided into two components: the quasi-static part, which corresponds to the weight that each contact point bears under static conditions; and the dynamic part, that results from the interaction between the two structures. In the railway case, despite the fact that in most cases the dynamic forces are less than 15% of the quasi-static component (Kruse and Popp, 2001; Katou et al., 2008), their consideration is of the greatest importance for the calculation of the response of the track-soil system, namely at remote positions (free-field). The calculation of the dynamic forces is addressed in the present section. As mentioned in the introduction to this chapter, the dynamic response of the vehicle is caused, among other aspects, by longitudinal variations of the stiffness of the supporting structure and by geometric irregularities at the contact between the vehicle and the structure. Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 137 Other aspects that influence the dynamic behavior of the vehicle-track system are the geometry and the dynamic properties of the track and vehicle, and the travel speed of the vehicle. All these aspects must be considered when developing a numerical model for the solution of this problem. In this work, since it is assumed that the track-soil system is invariant in the longitudinal direction, the variations of the stiffness cannot be considered and consequently the excitation due to the discrete sleeper support is not accounted for. It must be mentioned, however, that there are works in the literature in which the periodicity of the track-soil system is considered (Gupta et al., 2007; Gupta et al., 2008; Gupta et al., 2010). The geometric irregularities at the contact between the vehicle and the supporting structure (in the railway case, wheel-rail contact) are considered by means of a position dependent gap/irregularity profile. The irregularity profile may account for the unevenness of the track and for the imperfections of the wheels (Wu and Thompson, 2002). Besides the assumption of invariance in the longitudinal direction, this work also assumes that the vehicle and the supporting structure are linear and that the contact between the two structures is also linear. These two assumptions, which are commonly used by other authors (Metrikine et al., 2005; Lombaert et al., 2006), allow the analyses to be performed in the frequency domain, which is very convenient, since the 2.5D BEM-FEM yields results in the wavenumber-frequency domain and their transformation to the space-frequency domain is straightforward. Before proceeding to the description of the solution method, it must be mentioned that if the non-linear behavior of the vehicle or track is to be considered, or if the loss of contact between the wheels and the rail is to be studied, then frequency domain procedures lose their applicability and therefore time domain procedures must be used (Lane et al., 2007; Katou et al., 2008; Neves et al., 2012). Hybrid methods, in which the response of the track-soil system is first obtained in the 2.5D domain, then transformed to the space-time domain, and finally given as inputs to time domain procedures, can also be applied (Grundmann and Lenz, 2003; Müller et al., 2008). Nevertheless, the increase of complexity associated with these approaches is considerable and would force simplifications in other components of the system, namely in the boundary conditions of the track and length of the model. In the next sub-sections, the procedure used to calculate the dynamic forces that the vehicle transmits to the supporting structure is explained and exemplified. 4.3.2 Vehicle – structure interaction Consider a vehicle moving with constant speed V on top of an invariant structure, as represented in Figure 4.25a. The vehicle contacts with the structure through CP N contact points which may or may not belong to the same longitudinal alignment (in the example shown in Figure 4.25, the contact points belong to two distinct horizontal alignments). Chapter 4 – Invariant structures subjected to moving loads and moving vehicles 138 Figure 4.25: Vehicle-structure interaction – cross-section perspective (left) and longitudinal perspective (right): a) coupled system; b) substructuring method The solution method for the vehicle-structure interaction consists in a substructuring technique (Figure 4.25b), in which the vehicle and the supporting structure are modeled independently and in which the equilibrium of forces and the compatibility of displacements is enforced between the CP N contact points of the vehicle and the corresponding points of the supporting structure. The objective is to find the moving forces ( ) i f t ( CP 1... i N =) that the vehicle transmits to the supporting structure. Once the forces ( ) i f t are known, the response fields both of the vehicle and of the supporting structure can be calculated using the governing equations of each of the domains. The compatibility of displacements is imposed through the condition ( ) ( ) ( ) v s 0 0 , i i i i i u t u y Vt t u y Vt δ = + + + (4.32) in which ( ) v i u t is the displacement of the th i contact point of the vehicle, ( ) s 0 , i i u y Vt t + is the displacement of the corresponding moving contact point of the supporting structure and ( ) 0 i i u y Vt δ + is the irregularity/gap experienced by the contact point. The quantities s i u and i u δ are defined in a fixed frame of reference, and so 0 i y y Vt = + represents the longitudinal position of the th i contact point at the generic instant t , while 0 i y represents the corresponding position at the instant 0 t = . The displacements of the structure ( ) s 0 , i i u y Vt t + result from the contribution of all CP N contact points of the vehicle, and so they are calculated by the summation x z y z V a) b) ( ) i f t ( ) i f t V Analysis and mitigation of vibrations induced by the passage of high-speed trains in nearby buildings 139 ( ) ( ) ( ) ( ) CP 0 s s 0 0 1 , , j j N i i ij i f t y y Vt j u y Vt t u y Vt t δ − − + = + = + ∑ (4.33) where ( ) ( ) ( ) 0 s 0 , j j ij i f t y y Vt u y Vt t δ − − + + are the displacements that the force ( ) j f t − (force transmitted by the vehicle through the th j contact point) induces at the moving contact point i . The negative sign is used to account for the equilibrium condition. Since linear behavior and linear contact is assumed, eq. (4.32) can be transformed to the frequency domain, becoming ( ) ( ) ( ) v s 0 0 0 0 , i i i i u u y Vt u ω ω δ ω = + + ɶ ɶ ɶ (4.34) with ( ) ( ) 0 i v v 0 e d t i i u u t t ω ω +∞ − −∞ = ∫ ɶ (4.35) ( ) ( ) 0 i s s 0 0 0 , , e d t i i i i u y Vt u y Vt t t ω ω +∞ − −∞ + = + ∫ ɶ (4.36) ( ) ( ) 0 i 0 0 e d t i i i u u y Vt t ω δ ω δ +∞ − −∞ = + ∫ ɶ (4.37) Some comments must be made about these three frequency domain variables. The variable ( ) s 0 0 , i i u y Vt ω + ɶ represents the frequency domain displacements of a point with longitudinal coordinates 0 i y y Vt = + , i.e., a point that moves with the same speed as the set of loads. As mentioned in section 4.2.3 (equation (4.21)), the response field of points moving with the same speed as the load oscillate also with the same frequency as the load. Hence, accounting for (4.33) and (4.21), the integral (4.36) can be replaced with ( ) ( ) 0 0 00 0 CP ii s0 0 0 0 1 e , , e d 2 ij j i j i h y y y y NVV i i j ij j u y Vt f u V V ω ω ω ω ω ω ω ω π − −− +∞ =−∞ −   + = −     ∑∫ ɶ  ɶ ɶ ɶ (4.38) where ( ) 0 j f ω ɶ is the frequency content of the load transmitted by the th j contact point and calculated with ( ) ( ) 0 i 0 e d t j j f f t t ω ω +∞ − −∞ = ∫ ɶ (4.39) and where ( ) , ij y u k ω ɶ is the transfer function that relates the displacements of the alignment associated with the th i contact point with the forces at the alignment associated with the th j contact point. In this work, the referred to transfer function is calculated with the 2.5D BEMFEM procedure explained in chapter 3. On the other hand, the irregularity profile i u δ is a position dependent function, and so it is more convenient to define its Fourier transform in terms of the wavenumber y k rather than the radial frequency 0 ω . This leads to