# Numerical modeling of ground borne vibrations from underground traffic

Pranesh Chatterjee, Geert Degrande, Didier Clouteau, Tahmeed Al-Hussaini,
Maarten Arnst, Ramzi Othman

► To cite this version:

Pranesh Chatterjee, Geert Degrande, Didier Clouteau, Tahmeed Al-Hussaini, Maarten Arnst, et al. Numerical modeling of ground borne vibrations from underground traffic. 6th National Congress on theoretical and applied mechanics, 2003, Gent, Belgium. hal-01008546

HAL Id: hal-01008546

https://hal.science/hal-01008546

Submitted on 18 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Numerical modelling of ground borne vibrations from underground railway traffic

Pranesh Chatterjee¹, Geert Degrande¹, Didier Clouteau²,

Tahmeed Al-Hussaini², Maarten Arnst² and Ramzi Othman²

¹K.U.Leuven, Department of Civil Engineering, Structural Mechanics, Kasteelpark Arenberg 40, B-3001 Heverlee, Belgium

²Ecole Centrale de Paris, LMSSM, F-92295 Châtenay-Malabry, France
email: pranesh.chatterjee@bwk.kuleuven.ac.be, geert.degrande@bwk.kuleuven.ac.be,
clouteau@mss.ecp.fr, othman@mss.ecp.fr

*Abstract*— Vibrations due to the passage of trains in tunnels propagate through the soil and produce vibrations and re-radiated noise in adjacent structures. Within the frame of the EC-Growth project CONVURT ("The control of vibration from underground railway traffic"), an efficient and modular numerical prediction tool is being developed to predict vibration and re-radiated noise in adjacent buildings from excitation due to metro trains in tunnels for both new-build and existing situations. The model will be validated by means of in situ experiments at the site of Cité Universitaire in Paris and at Regent's Park in London. The development of a coupled periodic finite element-boundary element method and the results obtained from the in-situ experiments conducted in Paris are described in the present paper.

*Keywords*— soil structure interaction, finite element, boundary element, Floquet transform, Green-Floquet function

## I. INTRODUCTION

The passage of trains in tunnels generates ground-borne vibrations that propagate to adjacent structures and may produce vibrations and re-radiated noise. The amplitude of vibrations depends on several factors as the vehicle characteristics, the train speed, the track and wheel irregularities, the properties of the tunnel, the propagation of waves through the soil and the properties of the structures.

Two-dimensional finite element models with (local) absorbing boundary conditions are frequently used [1], [2], [3], but necessitate important simplifications to translate the three-dimensional moving load into an equivalent line load and do not allow to incorporate three-dimensional structures. Moreover, two-dimensional models underestimate radiation damping in the soil. The computational results may give a qualitative indication of vibration isolation efficiency (insertion loss) in the design phase. Underlying simplifying assumptions, however, compromise their use as absolute or quantitative predictions.

Within the frame of the EC-Growth project CONVURT ("The control of vibration from underground railway traffic"), an efficient and modular numerical prediction tool is being developed to predict vibration and re-radiated noise in adjacent buildings from excitation due to metro trains in tunnels for both new-build and existing situations. The model will be validated by means of elaborate in situ experiments at a site on the RER B line of RATP at Cité Universitaire in Paris and at another site on the Bakerloo line of London Underground in Regent's Park, London.

The three-dimensional (3-D) dynamic tunnel-soil interaction problem is solved by a coupled substructure approach [4] exploiting the invariance or the periodicity of

the tunnel in the longitudinal direction [5], [6], [7]. The paper describes a 3-D numerical model based on a finite element formulation for the tunnel structure with a modal reduction technique coupled with a boundary element formulation for a horizontally layered soil medium. The soil-structure interaction problem for the whole tunnel-soil system is restricted within the analysis of a generic cell by using Floquet transform. The impedance of the soil along the interface between the tunnel and the soil is computed with a boundary element formulation, using the Green-Floquet functions defined on a generic cell with the same width as the tunnel segment. Use is made of recent progress to model moving sources at the surface of a layered halfspace [8], [9]. The modal coordinates are calculated in the frequency-wavenumber domain. The dynamic representation theorem is used to compute the radiated wave field in the soil, while the inverse Floquet transformation is used to obtain the response in the frequency domain.

Elaborate in-situ experiments have been performed at the site of Cité Universitaire in Paris for the purpose of site characterization as well as validation of the numerical model. A spectral analysis of surface waves (SASW) test has been performed to obtain the thickness and the dynamic characteristics of the shallow layers of the soil. The vibrations at different points in the tunnel, in the free field and in a five-story reinforced concrete building have been measured during the passage of service trains, as well as a test train at variable speed. Furthermore, transfer functions between the track and the free field and the building have been recorded, using an impact hammer.

The main objective of the paper is to describe the development of the numerical prediction model and to comment on the experimental results obtained at the site of Cité Universitaire in Paris.

## II. NUMERICAL PREDICTION MODEL

In this section, the stepwise development of the numerical prediction model for a 3-D dynamic coupled linear tunnel-soil interaction problem is discussed. The problem geometry, the Floquet transformation using the periodicity of the problem in longitudinal direction, the basic equations of the generic cell and the solution of the generic problem by periodic coupled FEM-BEM approach are discussed in the following subsections.

### A. Problem geometry

A long tunnel is considered as an unbounded open set $Ω_t$, with given elastic properties and assumed to be embedded
in a layered elastic halfspace $\Omega$. The horizontally layered soil domain is denoted by $\Omega_s$, while $\Omega (=\Omega_t \cup \Omega_s)$ represents the global domain.

The tunnel $\Omega_t$ is assumed to be periodic in the longitudinal $y$-direction, defined by the unit vector $\mathbf{e}_y$. Hence, the generic cell in the tunnel may be defined by:

$$ \tilde{\Omega} = \{\mathbf{x} \in \Omega \mid 0 < \mathbf{x} \cdot \mathbf{e}_y < L\} \quad (1) $$

The tunnel $\Omega_t$ is constructed as an assemblage of an infinite number of elementary bounded cells $\tilde{\Omega}_{tn}$. The definitions of $\Omega_t$ and $\tilde{\Omega}_{tn}$ are given as

$$ \Omega_t = \bigcup_{n=-\infty}^{+\infty} \tilde{\Omega}_{tn} \quad (2) $$

$$ \tilde{\Omega}_{tn} = \left\{ \mathbf{x} \in \Omega / \mathbf{x} - nL\mathbf{e}_y \in \tilde{\Omega}_t \right\} \quad (3) $$

The elementary soil domain $\tilde{\Omega}_s$ in the reference cell has a similar definition.

The interface between the tunnel $\Omega_t$ and the soil $\Omega_s$ is denoted by $\Sigma$. In the reference cell, the interface between the tunnel and the soil segments $\tilde{\Omega}_t$ and $\tilde{\Omega}_s$ is denoted by $\tilde{\Sigma}$. The total interface $\Sigma$ can be decomposed as follows:

$$ \Sigma = \bigcup_{n=-\infty}^{+\infty} \tilde{\Omega}_n \quad (4) $$

The boundary $\partial\tilde{\Omega}$ of the generic cell $\tilde{\Omega}$ can be decomposed as follows:

$$ \partial\tilde{\Omega} = \tilde{\Gamma}_{\sigma} \cup \tilde{\Gamma}_{\infty} \cup \Sigma_0 \cup \Sigma_L \quad (5) $$

where $\tilde{\Gamma}_{\sigma}$ is the part of the boundary where Neumann boundary conditions are specified. The left and the right boundaries of the reference cell are denoted by $\Sigma_0$ and $\Sigma_L$ and defined as follows:

$$ \begin{aligned} \Sigma_0 &= \left\{ \mathbf{x} \in \delta\tilde{\Omega} / \mathbf{x} \cdot \mathbf{e}_y = 0 \right\} \\ \Sigma_L &= \left\{ \mathbf{x} \in \delta\tilde{\Omega} / \mathbf{x} \cdot \mathbf{e}_y = L \right\} \end{aligned} \quad (6) $$

Sommerfeld's radiation conditions should hold on $\tilde{\Gamma}_{\infty}$.

## B. Floquet transformations

The problem of dynamic soil-tunnel interaction can be considered as a three-dimensional problem having a periodicity with period $L$ in one direction $\mathbf{e}_y$. The position vector $\mathbf{x} \in \mathbb{R}^3$ of any point in the problem domain $\Omega$ is conveniently decomposed as:

$$ \mathbf{x} = \mathbf{x}' + ye_y = \mathbf{x}' + \tilde{y}\mathbf{e}_y + nL\mathbf{e}_y = \tilde{\mathbf{x}} + nL\mathbf{e}_y \quad (7) $$

where $\mathbf{x}'$ denotes the position vector in the $(x, z)$-plane perpendicular to the longitudinal direction $\mathbf{e}_y$, while $\tilde{\mathbf{x}}$ is the position vector in the reference cell $\tilde{\Omega}$.

The Floquet transform $\tilde{f}(\tilde{\mathbf{x}}, \kappa)$ of a general non-periodic function $f(\mathbf{x}) = f(\tilde{\mathbf{x}} + nL\mathbf{e}_y)$ defined on a three-dimensional problem domain $\tilde{\Omega}$ which is periodic with period $L$ in the direction $\mathbf{e}_y$ is defined as:

$$ \tilde{f}(\tilde{\mathbf{x}}, \kappa) = \sum_{n=-\infty}^{+\infty} f(\tilde{\mathbf{x}} + nL\mathbf{e}_y) \exp(+inL\kappa) \quad (8) $$

$\tilde{f}(\tilde{\mathbf{x}}, \kappa)$ is a continuous function of $\kappa$ in the interval $\kappa \in [-\pi/L, +\pi/L]$. It is a periodic function of the second kind with respect to $\tilde{\mathbf{x}}$ or $\tilde{y}$ as the following relation holds:

$$ \tilde{f}(\mathbf{x}', L, \kappa) = \exp(-i\kappa L)\tilde{f}(\mathbf{x}', 0, \kappa) \quad (9) $$

If the Floquet transform $\tilde{f}(\tilde{\mathbf{x}}, \kappa)$ is known on the reference cell $\tilde{\Omega}$, the function $f(\tilde{\mathbf{x}} + nL\mathbf{e}_y)$ can be reconstructed for any $\mathbf{x} = \tilde{\mathbf{x}} + nL\mathbf{e}_y$ using the inverse Floquet transform:

$$ f(\tilde{\mathbf{x}} + nL\mathbf{e}_y) = \frac{L}{2\pi} \int_{-\pi/L}^{+\pi/L} \tilde{f}(\tilde{\mathbf{x}}, \kappa) \exp(-inL\kappa)d\kappa \quad (10) $$

Equations (8) and (10) define the Floquet transform pair.

## C. Generic problem

It has been discussed in the previous subsection that, for a three-dimensional problem defined on a three-dimensional domain with periodicity $L$ along longitudinal direction $\mathbf{e}_y$, the analysis can be restricted to the generic domain $\tilde{\Omega}$. The boundary of $\partial\tilde{\Omega}$ of the generic domain $\tilde{\Omega}$ is decomposed as $\partial\tilde{\Omega} = \tilde{\Gamma}_{\sigma} \cup \Sigma_0 \cup \Sigma_L$ where $\tilde{\Gamma}_{\sigma}$ is the restriction of $\Gamma_{\sigma}$ on the generic cell and $\Sigma_0$ and $\Sigma_L$ are additional boundaries on which periodic conditions are imposed. Specifying the equilibrium equations in $\tilde{\Omega}_s$ and $\tilde{\Omega}_t$, the boundary conditions on $\tilde{\Gamma}_{\sigma s}$ and $\tilde{\Gamma}_{\sigma t}$ and the coupling equations on $\tilde{\Sigma}$, the following generic problem have to be solved.

For every $\kappa \in [-\pi/L, +\pi/L]$ and for every frequency $\omega \in \mathbb{R}$, the soil displacements $\tilde{\mathbf{u}}_s(\tilde{\mathbf{x}}, \kappa, \omega)$ in $\tilde{\Omega}_s$ should satisfy the following equations:

$$ \operatorname{div} \sigma_s(\tilde{\mathbf{u}}_s) = -\rho^s \omega^2 \tilde{\mathbf{u}}_s \quad \text{in} \quad \tilde{\Omega}_s \quad (11) $$

$$ \mathbf{t}_s(\tilde{\mathbf{u}}_s) = \mathbf{0} \quad \text{on} \quad \tilde{\Gamma}_{\sigma s} \quad (12) $$

$$ \tilde{\mathbf{u}}_s(\tilde{\mathbf{x}}) = \exp(-i\kappa L)\tilde{\mathbf{u}}_s(\tilde{\mathbf{x}} - L\mathbf{e}_y) \quad \text{on} \quad \tilde{\Sigma}_{Ls} \quad (13) $$

The tunnel displacements $\tilde{\mathbf{u}}_t(\tilde{\mathbf{x}}, \kappa, \omega)$ in $\tilde{\Omega}_t$ should satisfy the following equations:

$$ \operatorname{div} \sigma_t(\tilde{\mathbf{u}}_t) + \rho^t \mathbf{b} = -\rho^t \omega^2 \tilde{\mathbf{u}}_t \quad \text{in} \quad \tilde{\Omega}_t \quad (14) $$

$$ \mathbf{t}_t(\tilde{\mathbf{u}}_t) = \mathbf{f}_t \quad \text{on} \quad \tilde{\Gamma}_{\sigma t} \quad (15) $$

$$ \tilde{\mathbf{u}}_t(\tilde{\mathbf{x}}) = \exp(-i\kappa L)\tilde{\mathbf{u}}_t(\tilde{\mathbf{x}} - L\mathbf{e}_y) \quad \text{on} \quad \tilde{\Sigma}_{Lt} \quad (16) $$

On the boundary $\tilde{\Sigma}$ between the tunnel $\tilde{\Omega}_t$ and the soil $\tilde{\Omega}_s$, displacement continuity and stress equilibrium conditions are imposed by the following equations:

$$ \tilde{\mathbf{u}}_s = \tilde{\mathbf{u}}_t \quad \text{on} \quad \tilde{\Sigma} \quad (17) $$

$$ \mathbf{t}_s(\tilde{\mathbf{u}}_s) + \mathbf{t}_t(\tilde{\mathbf{u}}_t) = \mathbf{0} \quad \text{on} \quad \tilde{\Sigma} \quad (18) $$

The global solution is recovered using the inverse Floquet transform.

## D. Numerical solution of the generic problem

The numerical solution of the problem stated in the previous subsection will be built using the classical domain decomposition approach using FEM for the structure and BEM for the soil.
## D.1 The subdomain approach

As the tunnel domain $\tilde{\Omega}_t$ is bounded, the displacement field $\tilde{u}_t(\tilde{x}, \kappa, \omega)$ in the tunnel can be decomposed on a basis of functions $\Psi_n(\tilde{x}, \kappa)(n = 1, \dots, N)$ that are periodic of the second kind:

$$ \tilde{u}_t(\tilde{x}, \kappa, \omega) = \sum_{n=1}^{N} \Psi_n(\tilde{x}, \kappa) c_n(\kappa, \omega) = \Psi(\tilde{x}, \kappa) \mathbf{c}(\kappa, \omega) \quad (19) $$

The soil displacements $\tilde{u}_s(\tilde{x}, \kappa, \omega)$ are subsequently written as the sum of an elastodynamic wave field $\tilde{u}_0(\tilde{x}, \kappa, \omega)$ that is equal to zero on the interface $\tilde{\Sigma}$ and a scattered wave field $\tilde{u}_{sc}(\tilde{x}, \kappa, \omega)$:

$$ \tilde{u}_s(\tilde{x}, \kappa, \omega) = \tilde{u}_0(\tilde{x}, \kappa, \omega) + \tilde{u}_{sc}(\tilde{x}, \kappa, \omega) \quad (20) $$

The displacement $\tilde{u}_0(\tilde{x}, \kappa, \omega)$ is written as the sum of the incident wave field $\tilde{u}_i(\tilde{x}, \kappa, \omega)$ and the locally diffracted wave field $\tilde{u}_{d0}(\tilde{x}, \kappa, \omega)$, so that equation (20) may be rewritten as:

$$ \tilde{u}_s(\tilde{x}, \kappa, \omega) = \tilde{u}_i(\tilde{x}, \kappa, \omega) + \tilde{u}_{d0}(\tilde{x}, \kappa, \omega) + \tilde{u}_{sc}(\tilde{x}, \kappa, \omega) \quad (21) $$

The sum of the locally diffracted wave field and the scattered wave field is called the diffracted wave field. The displacement field $\tilde{u}_s(\tilde{x}, \kappa, \omega)$ in the soil can be further decomposed into:

$$ \tilde{u}_s(\tilde{x}, \kappa, \omega) = \tilde{u}_i(\tilde{x}, \kappa, \omega) + \tilde{u}_{d0}(\tilde{x}, \kappa, \omega) + \sum_{n=1}^{N} \tilde{u}_{dn}(\tilde{x}, \kappa, \omega)c_n(\kappa, \omega) \quad (22) $$

The following boundary conditions hold on the restriction $\tilde{\Sigma}$ of the interface $\Sigma$ on the reference cell:

$$ \begin{array}{ll} \tilde{u}_{dn}(\tilde{x}, \kappa) = \Psi_n(\tilde{x}, \kappa) & \text{on } \tilde{\Sigma} \\ \tilde{u}_i(\tilde{x}, \kappa, \omega) + \tilde{u}_{d0}(\tilde{x}, \kappa, \omega) = \mathbf{0} & \text{on } \tilde{\Sigma} \end{array} \quad (23) $$

A standard Galerkin procedure is used to write the stress equilibrium on $\tilde{\Sigma}$ in a weak sense for any function $\Psi_n(\tilde{x}, \kappa)$ in the basis and the following linear system of equations is obtained:

$$ \begin{aligned} & [\mathbf{K}_t(\kappa) - \omega^2 \mathbf{M}_t(\kappa) + \mathbf{K}_s(\kappa, \omega)] \mathbf{c}(\kappa, \omega) \\ & = \mathbf{F}_t(\kappa, \omega) + \mathbf{F}_s(\kappa, \omega) \end{aligned} \quad (24) $$

where the matrices $\mathbf{K}_t(\kappa)$ and $\mathbf{M}_t(\kappa)$ are the stiffness and mass matrices of the tunnel, the matrix $\mathbf{K}_s(\kappa, \omega)$ is the impedance matrix of the soil, the vector $\mathbf{F}_t(\kappa, \omega)$ is the generalized force vector transmitted from the track to the tunnel and the vector $\mathbf{F}_s(\kappa, \omega)$ is the generalized force vector due to the incident wave field.

## D.2 The tunnel impedance

The kinematics of the structure is described by the dynamic eigenmodes $\Phi_n(\tilde{x})$ of the unit cell, that are the solutions of the generalized eigenvalue problem:

$$ \mathbf{K}_t^{\text{FEM}} \Phi_n = \omega_n^2 \mathbf{M}_t^{\text{FEM}} \Phi_n \quad (25) $$

where $\mathbf{K}_t^{\text{FEM}}$ and $\mathbf{M}_t^{\text{FEM}}$ are the finite element stiffness and mass matrix of the unit cell. The eigenmodes $\Phi_n(\tilde{x})$ are periodic functions of the first kind on the unit cell $\tilde{\Omega}_t$ as the following kinematical conditions are imposed:

$$ \Phi_n(\tilde{x} + L\mathbf{e}_y) = \Phi_n(\tilde{x}) \quad \forall \tilde{x} \in \tilde{\Sigma}_{Lt} \ \& \ \forall \tilde{x} + L\mathbf{e}_y \in \tilde{\Sigma}_{Rt} \quad (26) $$

The fields $\Psi(\tilde{x}, \kappa)$ are subsequently constructed as the following periodic functions of the second kind:

$$ \Psi_n(\tilde{x}, \kappa) = \exp(-i\kappa\tilde{y})\Phi_n(\tilde{x}) \quad (27) $$

The matrices $\mathbf{K}_t(\kappa)$ and $\mathbf{M}_t(\kappa)$ in the governing equation (24) are the stiffness and mass matrices of the tunnel:

$$ K_{tIJ} = \int_{\tilde{\Omega}_t} \sigma(\Psi_I^*) : \epsilon(\Psi_J)d\Omega = \Psi_I^{*\mathrm{T}} \mathbf{K}_t^{\mathrm{FEM}} \Psi_J \quad (28) $$

$$ M_{tIJ} = \int_{\tilde{\Omega}_t} \rho_t \Psi_I^* \cdot \Psi_J d\Omega = \Psi_I^{*\mathrm{T}} \mathbf{M}_t^{\mathrm{FEM}} \Psi_J \quad (29) $$

while $\mathbf{F}_t$ is the generalized force vector transmitted from the track to the tunnel:

$$ F_{tI} = \int_{\tilde{\Sigma}} \Psi_I^* \cdot \tilde{f} dV = \Psi_I^{*\mathrm{T}} \tilde{f} \quad (30) $$

## D.3 The soil impedance

The matrix $\mathbf{K}_t(\kappa, \omega)$ in equation (24) is the impedance matrix of the soil, defined as:

$$ K_{sIJ} = \int_{\tilde{\Sigma}} \Psi_I^* \cdot \mathbf{t}_s(\tilde{u}_{dJ}) d\Sigma \quad (31) $$

where the stresses $\mathbf{t}_s(\tilde{u}_{dJ})$ are calculated with a periodic boundary element formulation with Green-Floquet functions defined on the periodic structure with period $L$ along the tunnel. In equation (24), $\mathbf{F}_s(\kappa, \omega)$ is the generalized force vector due to an incident wave field:

$$ F_{sI} = - \int_{\tilde{\Sigma}} \Psi_I^* \cdot (\mathbf{t}(\tilde{u}_i) + \mathbf{t}(\tilde{u}_{d0})) d\Sigma \quad (32) $$

The stress vector $\mathbf{t}(\tilde{u}_{d0})$ is also computed using a periodic boundary element formulation. In the case of dynamic tunnel-soil interaction problem where the external loading is applied on the tunnel invert, however, the force vector $\mathbf{F}_s(\kappa, \omega)$ can generally be disregarded.

If $u_{ij}^G(\mathbf{y}, \mathbf{x}, \omega)$ denotes the Green's tensor or the response in the direction $\mathbf{e}_j$ in the receiver $\mathbf{x}$ of a horizontally layered halfspace due to a unit harmonic force in the direction $\mathbf{e}_i$ in the source point $\mathbf{y}$, then the Green-Floquet tensor $\tilde{u}_{ij}^{GF}(\tilde{\mathbf{y}}, \tilde{\mathbf{x}}, \kappa, \omega)$ is defined as:

$$ \tilde{u}_{ij}^{GF}(\tilde{\mathbf{y}}, \tilde{\mathbf{x}}, \kappa, \omega) = \sum_{n=-\infty}^{+\infty} \exp(+in\kappa L) u_{ij}^{G}(\mathbf{y} + nL\mathbf{e}_y, \mathbf{x}, \omega) \quad (33) $$

The dynamic representation theorem now becomes:

$$ \begin{aligned} & c(\tilde{\xi}) \tilde{u}_i(\tilde{\xi}, \kappa, \omega) \\ & = \int_{\tilde{\Sigma} \cup \tilde{\Sigma}_{L_s} \cup \tilde{\Sigma}_{R_s}} \tilde{t}_j(\tilde{\mathbf{u}})(\tilde{x}, \kappa, \omega) \tilde{u}_{ij}^{GF}(\tilde{\xi}, \tilde{x}, \kappa, \omega) d\Sigma \\ & - \int_{\tilde{\Sigma} \cup \tilde{\Sigma}_{L_s} \cup \tilde{\Sigma}_{R_s}} \tilde{t}_{ij}^{GF}(\tilde{\xi}, \tilde{x}, \kappa, \omega) \tilde{u}_j(\tilde{x}, \kappa, \omega) d\Sigma \end{aligned} \quad (34) $$
After discretization of the unknown displacement and traction vectors along the boundary $\tilde{\Sigma}$ between the tunnel and the soil, the following system of equations is obtained:

$$ \tilde{\mathbf{U}}^{\mathrm{GF}}(\kappa, \omega) \tilde{\mathbf{t}}(\kappa, \omega) = \tilde{\mathbf{T}}^{\mathrm{GF}}(\kappa, \omega) \tilde{\mathbf{u}}(\kappa, \omega) \quad (35) $$

This system of equations is used to solve for the unknown tractions $\mathbf{t}_s(\tilde{\mathbf{u}}_{dJ})$ in equation (31), which allows to compute the soil impedance matrix $\mathbf{K}_t(\kappa, \omega)$.

### D.4 Wave propagation in the soil

Once the modal coordinates $\mathbf{c}(\kappa, \omega)$ are computed from the general equation (24), the displacements $\tilde{\mathbf{u}}_t(\tilde{\mathbf{x}}, \kappa, \omega)$ of the tunnel are obtained from equation (19). The tractions $\tilde{\mathbf{t}}_t(\Psi_n)(\tilde{\mathbf{x}}, \kappa, \omega)$ imposed by each mode $\Psi_n(\tilde{\mathbf{x}}, \kappa)$ are first computed using a periodic boundary element formulation. The traction vector $\tilde{\mathbf{t}}_t(\tilde{\mathbf{x}}, \kappa, \omega)$ on $\tilde{\Sigma}$ is finally obtained as:

$$ \tilde{\mathbf{t}}_t(\tilde{\mathbf{x}}, \kappa, \omega) = \sum_{n=1}^{N} \tilde{\mathbf{t}}_t(\Psi_n)(\tilde{\mathbf{x}}, \kappa, \omega) \tilde{c}_n(\kappa, \omega) \quad (36) $$

The incident wave field $\tilde{\mathbf{u}}^{\mathrm{inc}}(\tilde{\xi}, \kappa, \omega)$ is obtained by application of the dynamic representation theorem in the unbounded soil domain corresponding to the reference cell:

$$ \tilde{u}_i^{\mathrm{inc}}(\tilde{\xi}, \kappa, \omega) = \int_{\tilde{\Sigma}} \tilde{u}_{ij}^{\mathrm{GF}}(\tilde{\xi}, \tilde{\mathbf{x}}, \kappa, \omega) \tilde{\mathbf{t}}_{tj}(\tilde{\mathbf{x}}, \kappa, \omega) \\ - \tilde{\mathbf{t}}_{ij}^{\mathrm{GF}}(\tilde{\xi}, \tilde{\mathbf{x}}, \kappa, \omega) \tilde{\mathbf{u}}_{tj}(\tilde{\mathbf{x}}, \kappa, \omega) d\Sigma \quad (37) $$

where $u_{ij}^{\mathrm{GF}}(\tilde{\xi}, \tilde{\mathbf{x}}, \kappa, \omega)$ and $t_{ij}^{\mathrm{GF}}(\tilde{\xi}, \tilde{\mathbf{x}}, \kappa, \omega)$ are the Green-Floquet tensors and the integration is performed on the tunnel-soil interface $\tilde{\Sigma}$ in the reference cell.

The incident wave field in the soil is obtained after evaluation of the inverse Floquet transformation.

### III. NUMERICAL EXAMPLE

The metro tunnel on the line RER B of RATP at Cité Universitaire in Paris is a masonry cut-and-cover tunnel with two tracks at a shallow depth of about 9.3 m below the free surface and a width of 11.9 m (figure 1). The slab thickness is 0.6 m at the top and 0.4 m at the bottom, while the wall thickness is 1.5 m. The masonry is assumed to have a Young's modulus $E^t = 14000$ MPa, a Poisson's ratio $\nu^t = 0.15$, a density $\rho^t = 2400$ kg/m$^3$ and a hysteresis damping ratio $\beta^t = 0.02$.

The tunnel is considered to be infinitely long with a cross section that is invariant in the longitudinal $e_y$ direction. The length $L$ of the reference cell is chosen as 0.3 m. The reference cell is modelled with 8-node isoparametric brick elements (figure 2). Preliminary calculations only consider the transfer functions between the tunnel invert and the free field; therefore, the track is not yet included in the finite element model of the reference cell.

The eigenmodes $\Phi_n(\tilde{\mathbf{x}})$ of the reference cell of the tunnel are subjected to the periodicity condition (26). Figure 3 shows two modes of the reference cell that involve flexural motion in the transverse $(x, z)$-plane. The modes $\Psi(\tilde{\mathbf{x}}, \kappa)$ are subsequently constructed using the periodicity condition (27). As calculations will be made for excitation frequencies upto 80 Hz, 30 eigenmodes of the tunnel segment

<figure><img src="image_2.png" /><figcaption>Fig. 1. Cross section of the metro tunnel on the line RER B of RATP at Cité Universitaire.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>a. Cross section.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>b. Reference cell.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>a. Mode 6 at 14.0 Hz.</figcaption></figure>

<caption>b. Mode 8 at 36.9 Hz.</caption>

<caption>Fig. 3. Flexural modes of the reference cell of the tunnel.</caption>

(upto 225 Hz) are used to describe the kinematics of the tunnel.

The tunnel is embedded in a soil consisting of a shallow layer of 1.6 m of fill material, a layer of Beauchamp
sand with a thickness of approximately 3.2 m and a stiffer layer with a thickness of 7.8 m with marl and gravel on top of a halfspace consisting of chalk. A Spectral Analysis of Surface Waves (SASW) test has been performed to determine the thickness and the dynamic characteristics of the soil layers [10] and revealed the presence of a shallow layer with a thickness of $d = 1.4$ m and a shear wave velocity $C_s = 115$ m/s on top of a layer with $d = 3.0$ m and $C_s = 220$ m/s on a halfspace with $C_s = 315$ m/s. A Poisson's ratio $\nu^s = 0.4$, a density $\rho^s = 1700$ kg/m$^3$ and a material damping ratio $\beta^s = 0.05$ is assumed in all layers.

### A. Response in the wavenumber-frequency domain

The response of the tunnel-soil system due to a unit harmonic load on the tunnel invert in the frequency range between 0.1 Hz and 80.1 Hz are computed. The frequency step is equal to 2 Hz. Two forces are applied on both edges of the generic cell, resulting in a uniform wavenumber content. Computations are made in the wavenumber-frequency domain. For low frequencies upto 15 Hz, wavenumbers are sampled upto $\kappa^{\max} = 0.6$ m$^{-1}$ with a step $\Delta\kappa = 0.015$ m$^{-1}$, which is largely sufficient as the solution is shown to rapidly decay with increasing $\kappa$. For higher frequencies between 15 and 80 Hz, a sampling in the slowness domain corresponding to $p^{\max} = \kappa^{\max}/\omega = 0.004$ s/m and $\Delta p = 0.0001$ s/m is preferred.

<figure><img src="image_6.png" /><figcaption>Fig. 4. Real part of the diagonal elements $K_t(6, 6)$ and $K_t(8, 8)$ of the tunnel impedance matrix as a function of the slowness $p = \kappa/\omega$ and the frequency $\omega$.</figcaption></figure>

Figure 4 shows the real parts of the diagonal elements $K_t(6, 6)$ and $K_t(8, 8)$ of the tunnel impedance, corresponding to the flexible bending modes 6 and 8, as a function of the slowness $p = \kappa/\omega$ and the frequency $\omega$. At low values of the slowness and the frequency, the inertial term dominates the real part of the tunnel impedance, while at high values of the slowness and the frequency, the stiffness term dominates the tunnel impedance.

### B. Impedance of the soil

The impedance matrix $K_s(\kappa, \omega)$ of the soil is computed according to equation (31), using a periodic boundary element formulation. Figure 5 shows the real and the imaginary part of the elements $K_s(6, 6)$ and $K_s(8, 8)$ of the soil impedance matrix, corresponding to modes 6 and 8, as a function of the slowness $p = \kappa/\omega$ and the frequency $\omega$.

The (negative) imaginary part represents attenuation due to radiation and material damping in the soil. At zero slowness or infinite phase velocity in the longitudi-

<figure><img src="image_7.png" /><figcaption>Fig. 5. Real and imaginary parts of the diagonal elements $K_s(6, 6)$ and $K_s(8, 8)$ of the soil impedance matrix as a function of the slowness $p = \kappa/\omega$ and the frequency $\omega$.</figcaption></figure>

nal $y$-direction, the absolute value of the imaginary part increases, resulting in higher attenuation. For a fixed frequency, the absolute value of the imaginary part decreases for increasing slowness or decreasing phase velocity; it finally becomes very small, reflecting the absence of radiation damping in the soil at high values of the slowness.

### C. Response of the tunnel

The solution of the general system of equations (24) in the wavenumber-frequency domain finally results in the modal coordinates $c(\kappa, \omega)$. Figure 6 shows the modulus of the modal coordinates $c_6$ and $c_8$ as a function of the slowness $p = \kappa/\omega$ and the frequency $\omega$.

<figure><img src="image_8.png" /><figcaption>Fig. 6. Modulus of the modal coordinates $c_6$ and $c_8$ of the tunnel as a function of the slowness $p = \kappa/\omega$ and the frequency $\omega$.</figcaption></figure>

The response of the tunnel in the wavenumber-frequency domain is obtained through the synthesis of the contribution of all modes according to equation (19). The stresses along the tunnel-soil interface are evaluated and the dynamic representation theorem (34) is used to compute the incident wave field in the soil in the wavenumber-frequency domain.

### D. Response in the free field

In the present example, the response is computed in the soil at points on the free surface and at the level of the
tunnel base, at both sides of the tunnel upto a distance of 24 m from the central axis of the tunnel, as well as on a vertical line between the tunnel roof and the free surface.

Figure 7a shows the modulus of the transfer functions between the tunnel and 4 points at the free surface in the same cross section where the force is applied. A dominant peak around 14 Hz can clearly be observed. Figure 8 shows the real part of the displacements of the tunnel and the soil due to a harmonic loading on the tunnel invert at 14 Hz. It is clearly demonstrated that the soil above the tunnel and the tunnel roof move in phase. The peak at 14 Hz on figure 7a is due to resonance that can be explained with a single degree of freedom system where the tunnel acts as a spring and the soil above the tunnel represents the mass. This is confirmed in figure 7b where the displacements of the tunnel roof and the soil above the tunnel are shown to be equal.

<figure><img src="image_9.png" /><figcaption>Fig. 7. Modulus of the transfer functions for points (a) at the free surface and (b) at the tunnel floor, the tunnel roof and the free surface.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Real part of the displacements of the tunnel and the soil due to a harmonic excitation on the tunnel invert at 14 Hz at (a) $t = 0$, (b) $t = T/8$, (c) $t = 2T/8$, and (d) $t = 3T/8$.</figcaption></figure>

Figure 9 shows the real part of the displacements of the tunnel and the soil due to a harmonic loading on the tunnel invert at 80 Hz. Waves are propagating in the tunnel in the longitudinal and the transverse direction and radiated into the soil. The projection of the displacement field on the free surface in figure 10 clearly reveals elliptical wave fronts, as the waves that propagate in the transverse direction have

a lower velocity than the waves that propagate along the longitudinal direction. The latter can be demonstrated as guided waves inside the tunnel-soil system.

<figure><img src="image_11.png" /><figcaption>Fig. 9. Real part of the displacements of the tunnel and the soil due to a harmonic excitation on the tunnel invert at 80 Hz at (a) $t = 0$, (b) $t = T/8$, (c) $t = 2T/8$, and (d) $t = 3T/8$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Real part of the displacements at $t = 0$ at 80 Hz at the free surface.</figcaption></figure>

## IV. EXPERIMENTAL RESULTS

The results of the free field transfer function measurements between the track and the free field are presented here. Before going into the discussion about the results, the characteristics of the track and the experimental setup are briefly reviewed below.

### A. The characteristics of the track

Two classical ballast tracks supported on mono-block concrete sleepers are running in the tunnel. The mass per unit length of the UIC60 rail is 60 kg/m and the vertical bending stiffness is $6.4 \times 10^6$ Nm$^2$. The rail is supported on grooved rubber pads with a thickness of 9 mm which are resting on the sleepers. Mono-block concrete sleepers (VXP U61) are used. Each sleeper has a length of 2.27 m, a width of 0.29 m and a thickness of 0.14 m. The sleeper distance is 0.6 m. The mass of each sleeper is 200 kg and Young's modulus is 30 GPa. The thickness of the ballast layer is 0.4 m.

The direct rail receptance, measured with a weak static pre-load and a small dynamic force amplitude, shows a first resonance around 80-90 Hz, from which a ballast stiffness
per unit area of 100 MN/m³ and a loss factor of 0.80 can be estimated. The second resonance frequency around 600 Hz allows to derive a dynamic rail pad stiffness of about 375 MN/m and a loss factor of 0.50 (measured at 600 Hz and a small force amplitude, without a static pre-load). The pinned-pinned frequency of the track is about 1200 Hz, which is a common value for tracks with a UIC60 rail and a sleeper distance of 0.60 m. Measurements of direct sleeper receptance confirm the values of the natural frequencies of the track as determined by the direct rail receptance.

## B. The experimental setup

During the transfer function measurements, the excitation on the left rail of track 1 is generated with an impact hammer with a mass of 5.3 kg and a soft head as to excite in the low frequency range below 400 Hz. A load cell is placed on the impact hammer and measures the impact force on the rail head. Figure 11 shows the time history and frequency content of the vertical impact force. The peak value of the force is about 10.75 kN, while the frequency content is limited to the range below 400 Hz.

<figure><img src="image_13.png" /><figcaption>Fig. 11. (a) Time history and (b) frequency content of impact force.</figcaption></figure>

Figure 12 shows the position of accelerometers for measuring transfer functions between the track and the free field. The impact force and corresponding accelerations are measured in 16 channels for 14 impacts of the hammer.

## C. Response in the free field

Figure 13 shows the time history of the free field vertical accelerations as a function of the distance perpendicular to the tunnel for a hammer impact on the left rail of track 1. The function values in this seismogram have been rescaled so that the arrival times and the dispersive nature of the wave propagation can better be appreciated. On the contrary, the attenuation of the waves with increasing distance to the track is no longer visible.

Figures 14 and 15 show the time history and frequency content of the free field vertical acceleration on the top of the tunnel (0 m) and at 64 m from the tunnel. The maximum amplitude of the acceleration $a_z(0, 0, 0, t)$ in the near field (figure 14a) is about 0.017 m/s². The peak value of the acceleration $a_z(64, 0, 0, t)$ (figure 15a) in the far field is 0.001 m/s². These results demonstrate the decrease in the amplitude of the acceleration with increasing distance from the tunnel due to material and radiation damping in the soil. It is observed from figures 14b and 15b that the dominant frequency range is between 30 to 120 Hz. At

<figure><img src="image_14.png" /><figcaption>Fig. 12. Location of the measurement points.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13. Time history of the free field vertical accelerations as a function of the distance to the tunnel.</figcaption></figure>

a large distance from the tunnel, the high frequencies are attenuated due to material damping in the soil. A spurious frequency component at 110 Hz can be noticed.

<figure><img src="image_16.png" /><figcaption>Fig. 14. (a) Time history and (b) frequency content of the free field acceleration on the tunnel top (0 m).</figcaption></figure>

## D. Validation of the numerical model

Figure 16 shows the comparison of the computed and measured vertical seismograms in the free field. The computation has been performed with the numerical model described earlier in this paper. It is seen from figure 16 that at distances 4 m, 8 m, 12 m, 16 m, 24 m and 32 m, the computed vertical seismograms match reasonably well with experimentally obtained seismograms, exhibiting a similar
<figure><img src="image_17.png" /><figcaption>Fig. 15. (a) Time history and (b) frequency content of the free field acceleration at 64 m.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16. Comparison of computed and experimentally obtained vertical seismograms.</figcaption></figure>

dispersive nature of propagating waves in the soil medium and arrival times with increasing distance. Refracted longitudinal waves and Rayleigh waves are seen in the figure. For points farther away from the tunnel, the waves travel for longer distances in the third layer and then refract to the surface.

## V. CONCLUSION

The theoretical basis for the development of the numerical model for prediction of ground borne vibrations has been discussed in this paper. The model is based on a coupled periodic FE-BE method using a substructure approach, where the Floquet transformation is used to exploit the periodicity of the tunnel in longitudinal direction.

The experimental validation in this paper is limited to a preliminary comparison of the transfer functions between the tunnel and the free field and the computed response. A reasonable good correspondence of the arrival times and dispersive nature of the waves is observed.

## ACKNOWLEDGEMENTS

The results presented in this paper have been obtained within the frame of the EC-Growth project G3RD-CT-2000-00381 CONVURT ("The control of vibration from underground railway traffic"). The financial support of the European Community is gratefully acknowledged. In situ vibration measurements have been performed in collaboration with Vibratec and RATP.

## REFERENCES

[1] K.H. Chua, T. Balendra, and K.W. Lo, "Groundborne vibrations due to trains in tunnels," *Earthquake Engineering and Structural Dynamics*, vol. 21, no. 5, pp. 445-460, 1992.

[2] K.H. Chua, K.W. Lo, and T. Balendra, "Building response due to subway train traffic," *Journal of Geotechnical Engineering, Proceedings of the ASCE*, vol. 121, no. 11, pp. 747-754, 1995.

[3] W. Rücker and S. Said, "Einwirkung von U-Bahnerschütterungen auf Gebäude; Anregung, Ausbreitung und Abschirmung," in *Workshop Wave '94, Wave propagation and Reduction of Vibrations*, N. Chouw and G. Schmid, Eds., Rühr University, Bochum, Germany, December 1994, pp. 59-78.

[4] D. Aubry and D. Clouteau, "A subdomain approach to dynamic soil-structure interaction," in *Recent advances in Earthquake Engineering and Structural Dynamics*, V. Davidovici and R.W. Clough, Eds., pp. 251-272. Ouest Editions/AFPS, Nantes, 1992.

[5] H.E.M. Hunt, "Modelling of rail vehicles and track for calculation of ground vibration transmission into buildings," *Journal of Sound and Vibration*, vol. 193, no. 1, pp. 185-194, 1996.

[6] M.L. Elhabre, *Modélisation de l'interaction sismique sol-fluide-parois moulées suivant une approche périodique*, Ph.D. thesis, Laboratoire de Mécanique des Sols, Structures et Matériaux, Ecole Centrale de Paris, 2000.

[7] D. Clouteau, M.L. Elhabre, and D. Aubry, "Periodic BEM and FEM-BEM coupling," *Computational Mechanics*, vol. 25, pp. 567-577, 2000.

[8] F.C.P. de Barros and J.E. Luco, "Moving Green's functions for a layered visco-elastic half-space," Tech. Rep., Department of Applied Mechanics and Engineering Sciences of the University of California, La Jolla, 1992.

[9] D. Aubry, D. Clouteau, and G. Bonnet, "Modelling of wave propagation due to fixed or mobile dynamic sources," in *Workshop Wave '94, Wave propagation and Reduction of Vibrations*, N. Chouw and G. Schmid, Eds., Rühr University, Bochum, Germany, December 1994, pp. 109-121.

[10] L. Pyl and G. Degrande, "Determination of the dynamic soil characteristics with the SASW method at the site of Cité universitaire in Paris." Report BWM-2002-08, Department of Civil Engineering, Katholieke Universiteit Leuven, October 2002, CONVURT EC-Growth Project G3RD-CT-2000-00381.