# A numerical model for ground-borne vibrations from underground railway traffic based on a periodic finite element-boundary element formulation

Geert Degrande, Didier Clouteau, Ramzi Othman, Marteen Arnst, Hamdi Chebli, Ralf Klein, Pranesh K Chatterjee, B. Janssens

► To cite this version:

Geert Degrande, Didier Clouteau, Ramzi Othman, Marteen Arnst, Hamdi Chebli, et al.. A numerical model for ground-borne vibrations from underground railway traffic based on a periodic finite element-boundary element formulation. Journal of Sound and Vibration, 2006, 293 (3-5), pp.645-666. 10.1016/j.jsv.2005.12.023. hal-01004923

HAL Id: hal-01004923

https://hal.science/hal-01004923

Submitted on 10 Jun 2017

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
# A numerical model for ground-borne vibrations from underground railway traffic based on a periodic finite element-boundary element formulation

G. Degrandeª, D. Clouteauᵇ, R. Othmanᵇ, M. Arnstᵇ, H. Chebliᵇ, R. Kleinª,
P. Chatterjeeª, B. Janssensª

ªDepartment of Civil Engineering, KU Leuven, Kasteelpark Arenberg 40, B-3001 Leuven, Belgium

ᵇLMSSMat, Ecole Centrale de Paris, Grande Voie des Vignes, F-92295 Châtenay-Malabry, France

A numerical model is presented to predict vibrations in the free field from excitation due to metro trains in tunnels. The three-dimensional dynamic tunnel-soil interaction problem is solved with a subdomain formulation, using a finite element formulation for the tunnel and a boundary element method for the soil. The periodicity of the geometry in the longitudinal direction of the tunnel is exploited using the Floquet transform, limiting the discretization to a single-bounded reference cell. The responses of two different types of tunnel due to a harmonic load on the tunnel invert are compared, both in the frequency-wavenumber and spatial domains. The first tunnel is a shallow cut-and-cover masonry tunnel on the Paris metro network, embedded in layers of sand, while the second tunnel is a deep bored tunnel of London Underground, with a cast iron lining and embedded in the London clay.

## 1. Introduction

Underground trains create vibrations which are transmitted through the soil and interact with the foundations of adjacent buildings, resulting in disturbance from vibrations (1–80 Hz) and re-radiated noise (1–200 Hz).

Numerical models with varying degrees of sophistication are under development to predict vibrations from trains running in tunnels. This work has benefited to a large extent from recently developed models for vibrations due to trains running at grade on discretely supported and continuous tracks.

Dynamic axle loads are mainly caused by wheel and rail roughness, impact excitation due to discontinuities of the wheels and the rails and parametric excitation due to the periodic support of the rails on sleepers. State-of-the-art modelling of vibrations due to trains in tunnels assumes that the tunnel invert and the surrounding soil are relatively stiff with respect to the track, resulting in a two-step approach where the dynamic axle forces on the rails are derived from an uncoupled track model and subsequently applied on a track-tunnel-soil
model to compute the radiation of waves from the tunnel into the soil. The objective of the present paper is to
focus on the second problem.

Two-dimensional (plane strain) finite element (FE) models have been used to solve the dynamic (track-)tunnel-soil interaction model, where the infinite extension of the soil domain is represented by local absorbing boundary conditions [1-3]. Alternatively, two-dimensional coupled finite element-boundary element (FE-BE) formulations have been used [4,5], where the discretization effort is limited to the finite element mesh of the tunnel and the boundary element (BE) mesh, which can be restricted to the interface between the tunnel and the soil if the Green's functions of a layered half-space are used as an alternative to the full space solutions. However, two-dimensional models necessitate important simplifications to translate the three-dimensional (moving) load into an equivalent line load and underestimate radiation damping into the soil.

Three-dimensional FE models with local absorbing boundary conditions become prohibitively large for applications in the wide frequency range under consideration, as they require an extended finite element mesh together with a small element size. BE formulations partially solve this problem as they reduce the discretization effort. Stamos and Beskos [6,7] have used a three-dimensional BE formulation for the seismic analysis of large underground structures embedded in a homogeneous half-space. For applications to railway-induced vibrations in tunnels, however, involving the analysis at higher frequencies, three-dimensional coupled FE-BE formulations would also become very expensive from the computational point of view due to the large dimension of the fully populated BE matrices resulting from the coupling of all degrees of freedom on the interface between the tunnel and the soil. This has been the main motivation for many researchers to exploit the longitudinal invariance or periodicity of the track-tunnel-soil system for trains running in tunnels or the track-soil system for trains running at grade.

Forrest [8], Hussein and Hunt [9,10] and Hussein [11] have developed a semi-analytical insertion loss model, that accounts for the essential three-dimensional dynamic interaction between the track, the tunnel and the soil, as well as for wave propagation in the soil. The tunnel is modelled as a thin cylindrical shell surrounded by a homogeneous elastic medium. The model allows for the analysis of deep-bored tunnels embedded in a homogeneous soil. The equilibrium equations are formulated in the wavenumber-frequency domain, exploiting the invariance of the geometry along the longitudinal axis of the tunnel; models of this type are therefore referred to as two-and-a-half dimensional (2.5D), although it must be understood that they capture the three-dimensional behaviour of the invariant tunnel-track-soil system. Hussein [11] has also studied the vibration isolation efficiency of a floating slab in the tunnel, accounting for three different arrangements of the slab bearings.

For trains running at the surface, Sheng et al. [12,13], Madshus and Kaynia [14] and Lombaert et al. [15-17] have developed 2.5D models in the wavenumber-frequency domain that account for the invariance of the track in the longitudinal direction. These models are used to compute train-induced vibrations due to moving axle loads; the model of Lombaert et al. also accounts for the interaction between the train and the track.

In the context of the seismic analysis of sheet pile walls of infinite length, considerable progress has been made regarding the solution of dynamic soil-structure interaction problems defined on periodic domains, using coupled periodic FE-BE formulations [18-20]. Within the frame of the EC-Growth project CONVURT [21], this methodology has been used to develop a modular numerical model to predict vibrations from excitation due to metro trains in tunnels for both newly built and existing situations [22]. It is assumed that the computation of the dynamic axle loads on the rail can be decoupled from the solution of the dynamic track-tunnel-soil interaction problem. The present paper concentrates on the three-dimensional dynamic track-tunnel-soil interaction problem that is solved with a subdomain formulation, using a finite element formulation for the tunnel and a boundary element method for the soil. The periodicity of the tunnel and the soil is exploited using the Floquet transform, limiting the discretization to a single bounded reference cell of the tunnel. It has been demonstrated that the proposed approach can reuse the existing three-dimensional BE technology for layered media as the periodic Green's kernels have the same singularities as the three-dimensional Green's kernels [20]. This is a major advantage compared to the so-called 2.5D approaches, where a translation invariant model is assumed and for which the analysis of all singularities has to be repeated after the Fourier transform along the tunnel axis, except for surface tracks where the singular stress kernel vanishes [23].
The model will be validated by means of in situ experiments that have been performed at a site in the Cité Universitaire on the line Réseau Express Régional (RER B) of Régie Autonome des Transports Parisiens (RATP) in Paris [24–27] and a site in Regent's Park on the Bakerloo line of London Underground [28,29]. The tunnel in Paris is a shallow cut-and-cover masonry tunnel with two tracks, embedded in layers of sand, gravel and marl, while the tunnel in London is a deep cut-and-cover tunnel with a cast iron lining and a single track, embedded in a homogeneous layer of London clay.

After a brief review of the governing system of equations, details on the geometry and construction of both tunnels are presented in the present paper (the reader is referred to other papers [18,20,22] for a complete account of the formulation and solution of periodic dynamic soil-structure interaction problems). The response of both tunnels (and the surrounding soil) due to a harmonic load on the tunnel invert is considered, allowing to draw conclusions on the dynamic behaviour of both tunnel-soil systems. It is also demonstrated how the vibration isolation efficiency of a floating slab track can efficiently be computed using a Craig-Bampton substructuring technique [30]. Validation of the numerical model using the experimental results will be considered in a sequel publication.

The response for a harmonic load at different positions along the longitudinal axis in the reference cell (or the transfer functions) is needed to compute the response due to moving loads. This methodology has been outlined in general terms by Clouteau et al. [31], and has been applied to compute the response of both invariant [17] and periodic [32] surface tracks. The application of this methodology to the response of periodic tunnel-soil systems due to moving loads will also be discussed in a future publication.

## 2. Dynamic tunnel-soil interaction model

### 2.1. Problem outline

The analysis of the three-dimensional dynamic soil-tunnel interaction problem is restricted to straight tunnels that are periodic with a period $L$ in the longitudinal direction $e_y$ along the tunnel axis and embedded in a horizontally layered soil medium. The resulting problem is periodic and can be restricted to periodic fields of the second kind defined on a reference cell $\tilde{\Omega}$ (Fig. 1) by application of Floquet's theorem. The boundary $\partial\tilde{\Omega}$ of this domain is decomposed into the free surface $\tilde{\Gamma}_\sigma$ and the boundaries $\Sigma_0$ and $\Sigma_L$ on which periodic conditions are imposed. The generic cell $\tilde{\Omega}$ is decomposed into two subdomains: the soil $\tilde{\Omega}_s$ and the tunnel $\tilde{\Omega}_t$. The interface between these subdomains is denoted by $\tilde{\Sigma}_{ts}$. The boundary $\tilde{\Gamma}_{s\sigma}$ is the free surface of the soil, while a surface force $\tilde{f}_t$ is applied on $\tilde{\Gamma}_{ts}$ (Fig. 1).

The position vector $\mathbf{x}$ of any point in the problem domain $\Omega$ is decomposed as $\mathbf{x} = \tilde{\mathbf{x}} + nLe_y$, where $\tilde{\mathbf{x}}$ is the position vector in the reference cell $\tilde{\Omega}$ and $n$ is the cell number. The Floquet transformation $\tilde{f}(\tilde{\mathbf{x}}, \kappa)$ of a non-periodic function $f(\mathbf{x}) = f(\tilde{\mathbf{x}} + nLe_y)$ defined on a three-dimensional domain $\Omega$, that is periodic in the direction $e_y$ with period $L$, transforms the distance $nL$ between the $n$th cell and the reference cell $\tilde{\Omega}$ to the wavenumber

<figure><img src="image_2.png" /><figcaption>Fig. 1. Problem outline and notations.</figcaption></figure>
κ and is defined as [20]

$$ \tilde{f}(\tilde{\mathbf{x}}, \kappa) = \sum_{n=-\infty}^{+\infty} f(\tilde{\mathbf{x}} + n\mathbf{L}\mathbf{e}_y) \exp(+i\mathbf{n}\mathbf{L}\kappa). \quad (1) $$

The complex-valued function $\tilde{f}(\tilde{\mathbf{x}}, \kappa)$ is a periodic function of $\kappa$ with a period $2\pi/L$ (or a periodic function of the first kind with respect to $\kappa$), as the following periodicity condition holds:

$$ \tilde{f}\left(\tilde{\mathbf{x}}, \kappa + \frac{2\pi}{L}\right) = \tilde{f}(\tilde{\mathbf{x}}, \kappa). \quad (2) $$

This complex function $\tilde{f}(\tilde{\mathbf{x}}, \kappa)$ is also periodic of the second kind with a period $L$ in the direction $\mathbf{e}_y$, since the following condition holds for all $\tilde{\mathbf{x}}$:

$$ \tilde{f}(\tilde{\mathbf{x}} + \mathbf{L}\mathbf{e}_y, \kappa) = \exp(-i\kappa L)\tilde{f}(\tilde{\mathbf{x}}, \kappa). \quad (3) $$

The function $f(\tilde{\mathbf{x}} + n\mathbf{L}\mathbf{e}_y)$ can be reconstructed for any $\mathbf{x} = \tilde{\mathbf{x}} + n\mathbf{L}\mathbf{e}_y$ using the inverse Floquet transform:

$$ f(\tilde{\mathbf{x}} + n\mathbf{L}\mathbf{e}_y) = \frac{L}{2\pi} \int_{-\pi/L}^{+\pi/L} \tilde{f}(\tilde{\mathbf{x}}, \kappa) \exp(-i\mathbf{n}\mathbf{L}\kappa) d\kappa. \quad (4) $$

## 2.2. Navier equations

Using the Floquet transformation, all displacement and traction fields $\mathbf{u}(\mathbf{x}, \omega)$ and $\mathbf{t}(\mathbf{x}, \omega)$ defined on the periodic domain $\Omega$ are transformed to the fields $\tilde{\mathbf{u}}(\tilde{\mathbf{x}}, \kappa, \omega)$ and $\tilde{\mathbf{t}}(\tilde{\mathbf{x}}, \kappa, \omega)$ defined on the generic cell $\tilde{\Omega}$. Subscripts $s$ and $t$ will be used in the following to denote displacements, tractions and other field variables related to the soil and tunnel subdomain, respectively.

The Navier equations in the soil domain $\tilde{\Omega}_s$ and the boundary conditions on $\tilde{\Gamma}_{s\sigma}$ are written as follows for every frequency $\omega \in \mathbb{R}$ and wavenumber $\kappa \in [-\pi/L, +\pi/L]$:

$$ \mathrm{div} \, \tilde{\sigma}_s(\tilde{\mathbf{u}}_s) = -\rho^s \omega^2 \tilde{\mathbf{u}}_s \quad \text{in} \quad \tilde{\Omega}_s, \quad (5) $$

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

$$ \tilde{\mathbf{u}}_s(\tilde{\mathbf{x}}) = \exp(-i\kappa L)\tilde{\mathbf{u}}_s(\tilde{\mathbf{x}} - \mathbf{L}\mathbf{e}_y) \quad \text{on} \quad \Sigma_L, \quad (7) $$

where Eq. (5) expresses the conservation of momentum (with $\rho^s$ the density of the soil material) of the soil in the reference cell, Eq. (6) represents zero stresses on the part $\tilde{\Gamma}_{s\sigma}$ of the free surface and Eq. (7) imposes the periodicity condition of the second kind on the soil displacement vector $\tilde{\mathbf{u}}_s(\tilde{\mathbf{x}})$. This displacement field should also satisfy Sommerfeld's radiation conditions.

The following Navier equations and boundary conditions hold in the tunnel domain $\tilde{\Omega}_t$ with boundary $\tilde{\Gamma}_{t\sigma}$:

$$ \mathrm{div} \, \tilde{\sigma}_t(\tilde{\mathbf{u}}_t) = -\rho^t \omega^2 \tilde{\mathbf{u}}_t \quad \text{in} \quad \tilde{\Omega}_t, \quad (8) $$

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

$$ \tilde{\mathbf{u}}_t(\tilde{\mathbf{x}}) = \exp(-i\kappa L)\tilde{\mathbf{u}}_t(\tilde{\mathbf{x}} - \mathbf{L}\mathbf{e}_y) \quad \text{on} \quad \Sigma_L, \quad (10) $$

where $\rho^t$ is the density of the tunnel material and $\tilde{\mathbf{f}}_t$ is the force per unit area applied on the boundary $\tilde{\Gamma}_{t\sigma}$. Similar remarks can be made as for the Navier equations and boundary equations for the soil subdomain.

Continuity of displacements and equilibrium of stresses must hold on the tunnel-soil interface $\tilde{\Sigma}_{ts}$:

$$ \tilde{\mathbf{u}}_t = \tilde{\mathbf{u}}_s \quad \text{on} \quad \tilde{\Sigma}_{ts}, \quad (11) $$

$$ \tilde{\mathbf{t}}_t(\tilde{\mathbf{u}}_t) + \tilde{\mathbf{t}}_s(\tilde{\mathbf{u}}_s) = \mathbf{0} \quad \text{on} \quad \tilde{\Sigma}_{ts}. \quad (12) $$
## 2.3. Weak variational formulation

The Navier equations (9) on $\tilde{\Omega}_t$ and the boundary equations on the boundary $\tilde{\Gamma}_{t\sigma}$ are multiplied with the complex conjugate $\tilde{v}_t(\tilde{\mathbf{x}}, \kappa)$ of any virtual field $\tilde{v}_t(\tilde{\mathbf{x}}, \kappa)$ and integrated over the problem domain $\tilde{\Omega}_t$ and the boundary $\tilde{\Gamma}_{t\sigma}$. Integration by parts on the volume integral containing the divergence of the stress tensor allows to derive the following weak variational form for the tunnel [18,20,22]:

$$ \int_{\tilde{\Omega}_t} \tilde{\mathbf{\epsilon}}(\tilde{\mathbf{v}}_t) : \tilde{\mathbf{\sigma}}(\tilde{\mathbf{u}}_t) dV - \omega^2 \int_{\tilde{\Omega}_t} \rho^t \tilde{\mathbf{v}}_t \cdot \tilde{\mathbf{u}}_t dV = \int_{\partial \tilde{\Omega}_t} \tilde{\mathbf{v}}_t \cdot \tilde{\mathbf{t}}_t(\tilde{\mathbf{u}}_t) dS + \int_{\tilde{\Gamma}_{t\sigma}} \tilde{\mathbf{v}}_t \cdot \tilde{\mathbf{t}}_t dS, \quad (13) $$

where $\tilde{\mathbf{\epsilon}}(\tilde{\mathbf{v}}_t) : \tilde{\mathbf{\sigma}}(\tilde{\mathbf{u}}_t)$ denotes the double contraction of the complex conjugate of the virtual strain tensor $\tilde{\mathbf{\epsilon}}(\tilde{\mathbf{v}}_t)$ and the stress tensor $\tilde{\mathbf{\sigma}}(\tilde{\mathbf{u}}_t)$, which is equivalent to the trace of the tensor product $\tilde{\mathbf{\epsilon}}(\tilde{\mathbf{v}}_t)^{\mathrm{T}} \tilde{\mathbf{\sigma}}(\tilde{\mathbf{u}}_t)$. In Eq. (13), the boundary $\partial\tilde{\Omega}_t$ consists of the tunnel-soil interface $\tilde{\Sigma}_{ts}$ and the two boundaries $\Sigma_0$ and $\Sigma_L$ at the two edges $\tilde{y} = \pm L/2$ of the generic cell. As the actual and the virtual displacement fields are periodic of the second kind, the contribution of the integral on the sum of the boundaries $\Sigma_0$ and $\Sigma_L$ vanishes [18]. Accounting for the stress equilibrium (12) along the tunnel-soil interface $\tilde{\Sigma}_{ts}$, the weak variational equation becomes

$$ \int_{\tilde{\Omega}_t} \tilde{\mathbf{\epsilon}}(\tilde{\mathbf{v}}_t) : \tilde{\mathbf{\sigma}}(\tilde{\mathbf{u}}_t) dV - \omega^2 \int_{\tilde{\Omega}_t} \rho^t \tilde{\mathbf{v}}_t \cdot \tilde{\mathbf{u}}_t dV + \int_{\tilde{\Sigma}_{ts}} \tilde{\mathbf{v}}_t \cdot \tilde{\mathbf{t}}_s(\tilde{\mathbf{u}}_{sc}(\tilde{\mathbf{u}}_t)) dS = \int_{\tilde{\Gamma}_{t\sigma}} \tilde{\mathbf{v}}_t \cdot \tilde{\mathbf{t}}_t dS, \quad (14) $$

where $\tilde{\mathbf{u}}_{sc}(\tilde{\mathbf{u}}_t)$ denotes the wave field that is scattered by the tunnel into the soil and that obeys the displacement continuity relation (11) on the tunnel-soil interface $\tilde{\Sigma}_{ts}$.

## 2.4. Coupled periodic FE–BE formulation

As the tunnel is bounded, the displacement field $\tilde{\mathbf{u}}_t(\tilde{\mathbf{x}}, \kappa, \omega)$ can be decomposed as the following linear combination of functions $\tilde{\psi}_m(\tilde{\mathbf{x}}, \kappa)$ that are periodic of the second kind:

$$ \tilde{\mathbf{u}}_t(\tilde{\mathbf{x}}, \kappa, \omega) = \sum_{m=1}^{M} \tilde{\psi}_m(\tilde{\mathbf{x}}, \kappa) \alpha_m(\kappa, \omega) = \tilde{\Psi}_t \alpha_t, \quad (15) $$

where the modal coordinates $\alpha_m(\kappa, \omega)$ are collected in a vector $\alpha_t$. The kinematic basis for the tunnel is determined such that the modes $\tilde{\psi}_m(\tilde{\mathbf{x}}, \kappa)$ are periodic of the second kind:

$$ \tilde{\psi}_m(\tilde{\mathbf{x}} + Le_y, \kappa) = \exp(-i\kappa L)\tilde{\psi}_m(\tilde{\mathbf{x}}, \kappa), \quad \forall \tilde{\mathbf{x}} \in \tilde{\Omega}_t | \tilde{\mathbf{x}} \cdot \mathbf{e}_y = -\frac{L}{2}. \quad (16) $$

The modes $\tilde{\psi}_m(\tilde{\mathbf{x}}, \kappa)$ are constructed as

$$ \tilde{\psi}_m(\tilde{\mathbf{x}}, \kappa) = \exp(-i\kappa e_y \cdot \tilde{\mathbf{x}}) \tilde{\psi}_m^0(\tilde{\mathbf{x}}), \quad (17) $$

where the vectors $\tilde{\psi}_m^0(\tilde{\mathbf{x}})$ represent the eigenmodes of the reference cell, which are enforced to be periodic of the first kind (for $\kappa = 0$):

$$ \tilde{\psi}_m^0(\tilde{\mathbf{x}} + Le_y) = \tilde{\psi}_m^0(\tilde{\mathbf{x}}), \quad \forall \tilde{\mathbf{x}} \in \tilde{\Omega}_t | \tilde{\mathbf{x}} \cdot \mathbf{e}_y = -\frac{L}{2}. \quad (18) $$

The number $M$ of modes $\tilde{\psi}_m^0(\tilde{\mathbf{x}}, \kappa)$ in the expansion (15) of the tunnel displacement vector must be sufficiently high so as to incorporate the contribution of the relevant lowest modes of the tunnel, which depends on the frequency content and spatial distribution of the loading on the tunnel invert. This study is the subject of a convergence analysis.

The soil displacements $\tilde{\mathbf{u}}_s(\tilde{\mathbf{x}}, \kappa, \omega)$ can be written as the superposition of waves that are radiated by the tunnel into the soil:

$$ \tilde{\mathbf{u}}_s(\tilde{\mathbf{x}}, \kappa, \omega) = \tilde{\mathbf{u}}_{sc}(\tilde{\mathbf{u}}_t)(\tilde{\mathbf{x}}, \kappa, \omega) = \sum_{m=1}^{M} \tilde{\mathbf{u}}_{sc}(\tilde{\psi}_m)(\tilde{\mathbf{x}}, \kappa, \omega) \alpha_m(\kappa, \omega). \quad (19) $$

The numerical solution of the dynamic tunnel-soil interaction problem is obtained using the classical domain decomposition approach based on the finite element method for the structure and the BE method for the soil.
The displacements in the structure $\tilde{u}_t(\tilde{x}, \kappa, \omega)$ are interpolated as

$$ \tilde{u}_t \approx N_t \tilde{u}_t = N_t \tilde{\Psi}_t \alpha_t, \quad (20) $$

where $N_t$ is a matrix of global finite element shape functions (of dimension $3 \times N$ in a three-dimensional continuum formulation with $N$ degrees of freedom) and $\tilde{u}_t$ is the $N$-dimensional vector with nodal degrees of freedom. Global shape functions are preferred to local element-based shape functions so as to obtain a concise mathematical derivation. In practice, however, global finite element matrices are obtained by the assembly of element matrices, based on locally defined shape functions. Employing the same interpolation for the virtual displacements $\tilde{v}_t(\tilde{x}, \kappa, \omega)$, the following system of equations is finally obtained:

$$ [\mathbf{K}_t(\kappa) - \omega^2 \mathbf{M}_t(\kappa) + \mathbf{K}_s(\kappa, \omega)] \boldsymbol{\alpha}(\kappa, \omega) = \mathbf{F}_t(\kappa, \omega), \quad (21) $$

where $\mathbf{K}_t(\kappa)$ is the projection of the stiffness matrix $\mathbf{K}_t^{\text{FE}}$ of the tunnel reference cell on the tunnel modes (which are periodic of the second kind):

$$ \mathbf{K}_t(\kappa) = \overline{\tilde{\Psi}_t}^{\text{T}} \mathbf{K}_t^{\text{FE}} \underline{\tilde{\Psi}_t} = \overline{\tilde{\Psi}_t}^{\text{T}} \int_{\Omega_t} (\mathbf{L} \mathbf{N}_t)^{\text{T}} \mathbf{D}(\mathbf{L} \mathbf{N}_t) \, \mathrm{d}V \, \underline{\tilde{\Psi}_t}. \quad (22) $$

In a three-dimensional continuum formulation, $\mathbf{L}$ is a $3 \times 6$ matrix with differential operators that relates the strain vector to the displacement vector, while $\mathbf{D}$ is a $6 \times 6$ matrix with constitutive coefficients that relates the stress vector to the strain vector. The finite element matrix $\mathbf{K}_t^{\text{FE}}$ of the reference cell is an $N$-dimensional matrix, whereas the dimension of the stiffness matrix $\mathbf{K}_t(\kappa)$ is equal to the number of modes $M$ in the decomposition (15).

The matrix $\mathbf{M}_t(\kappa)$ is the projection of the mass matrix $\mathbf{M}_t^{\text{FE}}$ of the tunnel reference cell on the tunnel modes:

$$ \mathbf{M}_t(\kappa) = \overline{\tilde{\Psi}_t}^{\text{T}} \mathbf{M}_t^{\text{FE}} \underline{\tilde{\Psi}_t} = \overline{\tilde{\Psi}_t}^{\text{T}} \int_{\Omega_t} \mathbf{N}_t^{\text{T}} \rho^t \mathbf{N}_t \, \mathrm{d}V \, \underline{\tilde{\Psi}_t}. \quad (23) $$

$\mathbf{F}_t(\kappa, \omega)$ is the generalized force vector applied on the tunnel invert:

$$ \mathbf{F}_t(\kappa, \omega) = \overline{\tilde{\Psi}_t}^{\text{T}} \int_{\tilde{\Omega}_t} \mathbf{N}_t^{\text{T}} \tilde{\mathbf{f}}_t \, \mathrm{d}S \quad (24) $$

and $\mathbf{K}_s(\kappa, \omega)$ is the dynamic stiffness matrix of the soil:

$$ \mathbf{K}_s(\kappa, \omega) = \int_{\tilde{\Omega}_t} \overline{\tilde{\Psi}_t}^{\text{T}} \mathbf{N}_t^{\text{T}} \tilde{\mathbf{t}}_s(\tilde{\mathbf{u}}_{sc}(\mathbf{N}_t \tilde{\Psi}_t)) \, \mathrm{d}S. \quad (25) $$

The stresses $\tilde{\mathbf{t}}_s(\tilde{\mathbf{u}}_{sc}(\mathbf{N}_t \tilde{\Psi}_t))$ on the tunnel-soil interface are calculated with a periodic boundary element formulation with Green-Floquet functions or periodic Green functions defined on the periodic structure with period $L$ along the tunnel. The reader is referred to complementary literature for more details on the periodic boundary element formulation [18–20,22].

## 2.5. Craig–Bampton substructuring method

In order to analyse the performance of different track structures with respect to vibrations generated in the free field and in adjacent buildings, it is advantageous to differentiate between the degrees of freedom of the tunnel invert and the track. The motivation for this approach is to avoid the recomputation of the soil impedance, which is only governed by the kinematics of the tunnel invert, when alternative track configurations are considered in the tunnel. Therefore, the displacement vector $\tilde{\mathbf{u}}_t(\tilde{x}, \kappa, \omega)$ of the tunnel is discretized as follows, using a Craig–Bampton substructuring method [30]:

$$ \tilde{\mathbf{u}}_t = [\mathbf{N}_{t_r} \quad \mathbf{N}_{t_f}] \begin{Bmatrix} \tilde{\mathbf{u}}_{t_r} \\ \tilde{\mathbf{u}}_{t_f} \end{Bmatrix} = [\mathbf{N}_{t_r} \quad \mathbf{N}_{t_f}] \begin{bmatrix} \overline{\tilde{\Psi}}_{t_r} & \overline{\tilde{\Psi}}_{t_f}^S \\ 0 & \overline{\tilde{\Psi}}_{t_f} \end{bmatrix} \begin{Bmatrix} \alpha_{t_r} \\ \alpha_{t_f} \end{Bmatrix}, \quad (26) $$
where the subscripts $t_r$ and $t_t$ refer to the track and the tunnel invert, respectively. The modes in Eq. (26) are periodic of the second kind and are constructed as follows from periodic modes of the first kind:

$$ \begin{bmatrix} \underline{\tilde{\Psi}}_{t_r} & \underline{\tilde{\Psi}}_{t_r}^s \\ \underline{0} & \underline{\tilde{\Psi}}_{t_t} \end{bmatrix} = \begin{bmatrix} \Lambda_{t_r t_r} & \underline{0} \\ \underline{0} & \Lambda_{t_t t_r} \end{bmatrix} \begin{bmatrix} \underline{\tilde{\Psi}}_{t_r}^0 & \underline{\tilde{\Psi}}_{t_r}^{s0} \\ \underline{0} & \underline{\tilde{\Psi}}_{t_t}^0 \end{bmatrix}, \qquad (27) $$

where the diagonal matrices $\Lambda_{t_r t_r}$ and $\Lambda_{t_t t_t}$ are constructed according to Eq. (17). The modes $\underline{\tilde{\Psi}}_{t_r}^0$ are the eigenmodes of the track clamped at the tunnel invert. The modes $\underline{\tilde{\Psi}}_{t_t}^0$ are the modes of the free tunnel without track. The displacements $\underline{\tilde{\Psi}}_{t_r}^{s0}$ are the quasi-static transmission of the tunnel modes into the track, computed as

$$ \underline{\tilde{\Psi}}_{t_r}^{s0} = -(\mathbf{K}_{t_r t_r}^{\mathrm{FE}})^{-1} \mathbf{K}_{t_r t_t}^{\mathrm{FE}} \underline{\tilde{\Psi}}_{t_t}^{0}, \qquad (28) $$

where $\mathbf{K}_{t_r t_r}^{\mathrm{FE}}$ and $\mathbf{K}_{t_t t_t}^{\mathrm{FE}}$ are block submatrices of the finite element stiffness matrix $\mathbf{K}^{\mathrm{FE}}$ of the reference cell of the tunnel. The Craig-Bampton decomposition is frequently used for the seismic analysis of structures on flexible foundations [33]. It is equivalent to a classical modal decomposition (based on the modes of the complete track-tunnel structure) if all modes are included in both formulations. If a reduced modal basis is used, a convergence analysis is needed to determine the number of tunnel and track modes that need to be included in the Craig-Bampton decomposition (26), analogously to a convergence analysis for classical modal analysis.

Introducing the decomposition (26) in Eq. (14) results in

$$ \left\{ \begin{bmatrix} \underline{\tilde{\Psi}}_{t_r} & \underline{\tilde{\Psi}}_{t_r}^s \\ \underline{0} & \underline{\tilde{\Psi}}_{t_t} \end{bmatrix}^T \left( \begin{bmatrix} \mathbf{K}_{t_r t_r}^{\mathrm{FE}} & \mathbf{K}_{t_r t_t}^{\mathrm{FE}} \\ \mathbf{K}_{t_t t_r}^{\mathrm{FE}} & \mathbf{K}_{t_t t_t}^{\mathrm{FE}} \end{bmatrix} - \omega^2 \begin{bmatrix} \mathbf{M}_{t_r t_r}^{\mathrm{FE}} & \mathbf{M}_{t_r t_t}^{\mathrm{FE}} \\ \mathbf{M}_{t_t t_r}^{\mathrm{FE}} & \mathbf{M}_{t_t t_t}^{\mathrm{FE}} \end{bmatrix} \right) \begin{bmatrix} \underline{\tilde{\Psi}}_{t_r} & \underline{\tilde{\Psi}}_{t_r}^s \\ \underline{0} & \underline{\tilde{\Psi}}_{t_t} \end{bmatrix} + \begin{bmatrix} \underline{0} & \underline{0} \\ \underline{0} & \mathbf{K}_s \end{bmatrix} \right\} \begin{Bmatrix} \alpha_{t_r} \\ \alpha_{t_t} \end{Bmatrix} = \begin{Bmatrix} \mathbf{F}_{t_r} \\ \mathbf{0} \end{Bmatrix}. \qquad (29) $$

$\mathbf{K}_s(\kappa, \omega)$ is the dynamic stiffness matrix of the soil:

$$ \mathbf{K}_s(\kappa, \omega) = \int_{\tilde{\Sigma}_{is}} \bar{\underline{\tilde{\Psi}}}_{t_r}^T \mathbf{N}_{t_r}^{\mathrm{T}} \tilde{\mathbf{t}}_s (\tilde{\mathbf{u}}_{sc}(\mathbf{N}_{t_r} \underline{\tilde{\Psi}}_{t_r})) \mathrm{d}\mathbf{S}. \qquad (30) $$

$\mathbf{F}_{t_r}(\kappa, \omega)$ is the generalized force vector applied on the track:

$$ \mathbf{F}_{t_r}(\kappa, \omega) = \bar{\underline{\tilde{\Psi}}}_{t_r}^T \int_{\tilde{\Gamma}_{t\sigma}} \mathbf{N}_{t_r}^T \tilde{\mathbf{f}}_{t_r} \mathrm{d}\mathbf{S}. \qquad (31) $$

The soil impedance in Eq. (30) only depends on the tunnel modes and does not change when a calculation needs to be made for an alternative track structure in the tunnel, which contributes to the versatility of the proposed approach.

## 2.6. Wave propagation in the soil

When the displacements $\tilde{\mathbf{u}}(\tilde{\mathbf{x}}, \kappa, \omega)$ and the stresses $\tilde{\mathbf{t}}(\tilde{\mathbf{x}}, \kappa, \omega)$ on the tunnel-soil interface are known, the incident wave field $\tilde{\mathbf{u}}^{\mathrm{inc}}(\tilde{\mathbf{\xi}}, \kappa, \omega)$ is obtained by application of the dynamic representation theorem in the unbounded soil domain corresponding to the reference cell:

$$ \tilde{\mathbf{u}}_i^{\mathrm{inc}}(\tilde{\mathbf{\xi}}, \kappa, \omega) = \int_{\tilde{\Sigma}_{is}} \tilde{\mathbf{u}}_{ij}^{\mathrm{GF}}(\tilde{\mathbf{\xi}}, \tilde{\mathbf{x}}, \kappa, \omega) \tilde{\mathbf{l}}_{ij}(\tilde{\mathbf{x}}, \kappa, \omega) - \tilde{\mathbf{l}}_{ij}^{\mathrm{GF}}(\tilde{\mathbf{\xi}}, \tilde{\mathbf{x}}, \kappa, \omega) \tilde{\mathbf{u}}_{ij}(\tilde{\mathbf{x}}, \kappa, \omega) \mathrm{d}\mathbf{S} \qquad (32) $$

with $\tilde{\mathbf{u}}_{ij}^{\mathrm{GF}}(\tilde{\mathbf{\xi}}, \tilde{\mathbf{x}}, \kappa, \omega)$ and $\tilde{\mathbf{l}}_{ij}^{\mathrm{GF}}(\tilde{\mathbf{\xi}}, \tilde{\mathbf{x}}, \kappa, \omega)$ the Green-Floquet or periodic Green tensors [18,20]. The incident wave field in the soil is obtained by evaluating the inverse Floquet transform (4).
## 3. Numerical results for the Bakerloo line tunnel in London

### 3.1. Site characteristics

The tunnel on the Bakerloo line of London Underground in Regent's Park is a deep bored tunnel with a cast iron lining and a single track, embedded in London clay at a depth of about 28 m. The tunnel has an internal radius of 1.83 m and a wall thickness of 0.022 m. There are six longitudinal stiffeners and one circumferential stiffener at an interval of 0.508 m, resulting in a periodic structure (Fig. 2). The track is a non-ballasted concrete slab track with Bull head rail supported on hard Jarrah wooden sleepers via cast iron chairs. The sleeper distance is 0.9 m. Both ends of a sleeper are concreted into the invert and the space between the sleepers is filled with shingle. Resilience is mainly provided by the timber sleepers, as the rails are not supported by rail pads.

Geological maps show that the thickness of the London clay layer at the site is 40 m. GeoDelft has performed cone penetration tests (CPT) upto a depth of 21 m [34]. A shallow top layer with a thickness of 4–6 m has inclusions of sand and gravel and varying cone resistance. Bender element tests on undisturbed samples at several confining pressures result in an average shear wave velocity of 124 m/s and a longitudinal wave velocity of 1604 m/s [34]. A material damping ratio of 0.042 in the top layer and 0.039 in the second layer has been determined with free torsion pendulum tests [34]. A seismic cone penetration test (SCPT) upto a depth of 21 m confirmed the presence of a shallow stiff layer with a thickness of 4–6 m and a shear wave velocity $C_s$ = 325 m/s on top of a homogeneous half-space with $C_s$ = 220 m/s [34]. Spectral analysis of surface wave (SASW) tests confirmed the latter and revealed the presence of a homogeneous clay layer with a shear wave velocity between 200 and 260 m/s [35].

In the numerical calculations, the tunnel is assumed to be embedded in a layered soil consisting of a single shallow layer with a thickness of 5 m on top of a homogeneous half-space consisting of clay. The top layer has

<figure><img src="image_3.png" /><figcaption>Fig. 2. Cross-section of the metro tunnel on the Bakerloo line at Regent's Park.</figcaption></figure>
a shear wave velocity $C_s = 275$ m/s, a longitudinal wave velocity $C_p = 1964$ m/s, a density $\rho^s = 1980$ kg/m$^3$ and a material damping ratio $\beta^s = 0.042$. The underlying half-space has a shear wave velocity $C_s = 220$ m/s, a longitudinal wave velocity $C_p = 1571$ m/s, a density $\rho^s = 1980$ kg/m$^3$ and a material damping ratio $\beta^s = 0.039$.

## 3.2. Kinematics of the tunnel

The tunnel in London is periodic, as its lining is made of cast iron segments with circumferential stiffeners. The length $L$ of the cell is equal to the segment length of 0.508 m. The sleeper distance is equal to 0.9 m and introduces a second periodicity, which is not accounted for as the track is not included in the model.

A FE model of the reference cell is made using the structural dynamics toolbox (Sdt) in Matlab. The exterior shell is modelled with four-node quadrilateral shell elements based on a Q4WT formulation for the membrane behaviour and the MITC4 formulation to represent the plate bending behaviour. In the latter, a mixed interpolation of the transverse displacement, section rotations and transverse shear strains is used to avoid shear locking when applied to thin shells [36]. The longitudinal and circumferential stiffeners are modelled with standard two-node Euler-Bernoulli beam elements, that are rigidly connected to the exterior shell using stiff eccentric connections. Eight-node brick elements are used to model the concrete on the tunnel invert.

The kinematic basis for the tunnel consists of the modes $\tilde{\psi}_m(\tilde{x}, \kappa)$ that are derived from the eigenmodes $\psi_m^0(\tilde{x})$ of the tunnel cell with free boundary conditions on $\tilde{\Sigma}_{Ts}$ and periodicity conditions at both ends $\Sigma_0$ and $\Sigma_L$. Due to these constraints and the symmetry of the cell, displacements in the $y$-direction are decoupled from displacements in the $x$- and $z$-direction and only four rigid body modes are found. Fig. 3 shows the first two in-plane and the first out-of-plane flexible modes $\psi_m^0(\tilde{x})$ of the reference cell.

<figure><img src="image_4.png" /><figcaption>Fig. 3. The first (a,b) in-plane and (c) out-of-plane modes of the reference cell of the Bakerloo line tunnel. (a) Mode 5 at 20.94 Hz, (b) mode 6 at 24.32 Hz and (c) mode 11 at 153.53 Hz.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Modulus of the element Z<sub>t</sub>(6, 6) of the tunnel impedance matrix of the Bakerloo line tunnel (a) as a function of ω and p and (b) as a function of ω for constant p = 0 s/m (solid line), p = 0.3 × 10<sup>-3</sup> s/m (dashed line) and p = 0.6 × 10<sup>-3</sup> s/m (dotted line). Both figures are on a logarithmic scale.</figcaption></figure>

In the following, the response of the tunnel-soil system due to an eccentric harmonic vertical load on the tunnel invert is computed in the frequency range upto 80 Hz. Convergence analysis has shown that minimum 20 tunnel modes of the reference cell need to be included in the analysis.

3.3. Impedance of the tunnel

Fig. 4a shows the modulus of the sixth diagonal element Z<sub>t</sub>(6, 6) of the tunnel impedance matrix Z<sub>t</sub>(κ, ω) = K<sub>t</sub>(κ) − ω<sup>2</sup>M<sub>t</sub>(κ) as a function of the frequency ω and the slowness p = κ/ω, which is defined as the ratio of the wavenumber κ to the circular frequency ω. The element Z<sub>t</sub>(6, 6) corresponds to the sixth mode of the reference cell (Fig. 3b) involving in-plane deformation of the reference cell, equivalent to a mode of azimuthal order m = 2 in an isotropic shell (i.e. with 2 wavelengths along the circumferential direction). Fig. 4a clearly illustrates that the sixth eigenfrequency of the tunnel reference cell corresponds to the cut-on frequency of a dispersive mode that propagates in the longitudinal direction of the tunnel. Fig. 4b shows the variation of the element Z<sub>t</sub>(6, 6) of the tunnel impedance matrix with the frequency ω for a zero value of the slowness p = 0 s/m, corresponding to the two-dimensional plane strain case. At ω = 0 rad/s, the tunnel impedance is equal to the square of the sixth eigenfrequency of the reference cell. Due to the inertial term, it decreases with the square of the frequency ω and becomes zero at the sixth resonance frequency of the reference cell, which corresponds to the cut-on frequency of a dispersive mode for which the phase velocity is infinite at zero slowness or wavenumber in the undamped case. For non-zero values of p, the tunnel impedance increases with the square of the slowness p. Fig. 4b also shows the variation of the element Z<sub>t</sub>(6, 6) of the tunnel impedance matrix with the frequency ω for two non-zero values of the slowness p = 0.3 × 10<sup>-3</sup> and 0.6 × 10<sup>-3</sup> s/m.

The response of the tunnel is determined by the contribution of different tunnel modes. Fig. 5 shows, on a
logarithmic scale, the modulus of the vertical displacement of the tunnel invert under the load as a function of
the circular frequency ω and the slowness p for the free tunnel that is not supported by the soil. In all, 20
modes of the reference cell are included in the analysis. The response is governed by dispersive tunnel modes,
that start at the eigenfrequencies of the tunnel reference cell; the contribution of the modes 5 till 8 with a cut-
on frequency at 20.94, 24.32, 62.74 and 68.04 Hz can clearly be observed.

3.4. Impedance of the soil

Fig. 6 shows the real and imaginary parts of the element K_s(6, 6) of the impedance matrix **K_s**(κ, ω) of the soil, corresponding to the sixth mode of the **Bakerloo** line tunnel, as a function of the frequency ω and the slowness *p*. Note that the range of the slowness *p* extends to larger values than was the case in Fig. 4a. The real
<figure><img src="image_6.png" /><figcaption>Fig. 5. Modulus of the vertical displacement (dB ref. $10^{-6}$ m) of the invert of the Bakerloo line tunnel as a function of $\omega$ and $p$, accounting for 20 tunnel modes.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. (a) Real and (b) imaginary parts of the element $K_x(6, 6)$ of the soil impedance matrix of the Bakerloo line tunnel as a function of $\omega$ and $p$.</figcaption></figure>

part increases for increasing values of $p$. Due to the sign convention used in the kernels of the Fourier and Floquet transformations, the imaginary part is negative. High (negative) values of the imaginary part reflect dissipation of energy due to radiation damping in the soil. The low imaginary part of the soil impedance for high values of the slowness $p$ (or decreasing phase velocity along the longitudinal tunnel axis) reflects the absence of radiation damping in the soil.

Boström and Burden [37] have demonstrated that the response of an infinite elastic medium with a cylindrical cavity is governed by dispersive surface waves modes, corresponding to integer values of the azimuthal wavenumber $m$; these modes are confined to the vicinity of the cavity and have a cut-on frequency depending on the azimuthal wavenumber $m$, the characteristics of the soil domain and the radius of the tunnel. They propagate with a phase velocity that varies between the shear wave velocity (at the cut-on frequency) and the Rayleigh wave velocity (at limiting high frequency). As the Bakerloo line tunnel is embedded in a deep homogeneous layer of clay, the behaviour of an element of the soil impedance matrix is very similar to the response of a cylindrical cavity in an elastic medium, when a particular mode of this cavity is considered. The difference, however, is that the kinematics along the tunnel-soil interface are governed by the modes of the (orthotropic) reference cell of the tunnel that do not perfectly match the azimuthal modes of a cylindrical cavity or an isotropic shell. Based on the analytical solution of Boström and Burden, the cut-on frequency of
the mode corresponding to the azimuthal mode number $m = 2$ of a cylindrical cavity with a radius of 1.953 m, embedded in a homogeneous and undamped medium with the characteristics of the London clay, is equal to 65 Hz.

### 3.5. Response due to harmonic loading

The response of the Bakerloo line tunnel is considered for a harmonic loading on the tunnel invert in the frequency range between 1 and 80 Hz.

Fig. 7a shows the modulus of the vertical displacement of the tunnel invert as a function of the circular frequency $\omega$ and the slowness $p$ for the case of a free tunnel that is not supported by the soil. These results correspond to those presented in Fig. 5, but now a wider range of the slowness $p$ is considered to enable the comparison with results that follow in this subsection. The response is governed by the dispersive character of the tunnel modes at low values of the slowness, as explained before.

Fig. 7b considers the response of the layered soil medium with a cylindrical cavity, only accounting for the impedance of the soil and neglecting the dynamic impedance of the tunnel. The kinematics of the free surface along the cavity (i.e. the interface between the soil and the tunnel in the dynamic soil-tunnel interaction problem) is described by 20 modes of the tunnel's reference cell. The response is dominated by dispersive

<figure><img src="image_8.png" /><figcaption>Fig. 7. Modulus of the vertical displacement (dB ref. $10^{-6}$ m) of the invert of the Bakerloo line tunnel as a function of $\omega$ and $p$, accounting for 20 tunnel modes. Results are shown for the case of (a) the free tunnel, (b) the layered soil with a cylindrical cavity, and (c) the coupled soil-tunnel problem.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 8. Modulus of the vertical displacement (dB ref. $10^{-6}$ m) of the invert of the Bakerloo line tunnel as a function of $p$ at (a) 21 Hz, (b) 41 Hz, (c) 61 Hz, and (d) 81 Hz, accounting for 20 tunnel modes. Results are shown for the case of the free tunnel (dashed line), the layered soil with a cylindrical cavity (dash-dotted line), and the coupled soil-tunnel problem (solid line).</figcaption></figure>

modes that propagate with a phase velocity that varies between the shear wave velocity $C_s = 220$ m/s of the soil around the tunnel (at the cut-on frequency of each mode) and the Rayleigh wave velocity $C_r = 209.75$ m/s at limiting high frequency. The corresponding values of the slowness vary between $1/C_s = 4.545 \times 10^{-3}$ s/m and $1/C_r = 4.768 \times 10^{-3}$ s/m.

Fig. 7c finally shows the response of the coupled soil-tunnel system in the frequency-slowness domain. Furthermore, Fig. 8 shows sections of the vertical response of the tunnel invert as a function of the slowness $p$ for four values of the frequency (21, 41, 61 and 81 Hz). Results are shown for the three cases considered earlier: the free tunnel (Fig. 7a), the soil medium with a cylindrical cavity (Fig. 7b) and the coupled soil-tunnel system (Fig. 7c). Fig. 8b, for example, clearly illustrates that the response of the free tunnel at 41 Hz is dominated by two dispersive bending modes of the tunnel, while the response of the soil medium with a cylindrical cavity is maximum for values of the slowness between the reciprocals of the shear wave velocity and Rayleigh wave velocity of the soil around the tunnel. For low values of the slowness, the dynamic stiffness of the soil is larger than the dynamic stiffness of the tunnel, and the peaks in the response of the free tunnel are no longer apparent in the response of the coupled tunnel-soil system. At high values of the slowness, the dynamic stiffness of the tunnel is very high compared to the soil's stiffness, and the response of the coupled soil-tunnel system converges to the response of the free tunnel, while the effect of dispersive modes along the tunnel-soil interface is no longer observable.

The response in the spatial domain is obtained after evaluation of the inverse Floquet transformation. Figs. 9 and 10 show the response of the Bakerloo line tunnel for a harmonic excitation on the tunnel invert at 10 and 40 Hz, respectively. The tunnel's cross section is relatively stiff but the tunnel diameter is small, contributing to the flexibility of the tunnel in the longitudinal direction. As the Bakerloo line tunnel is embedded at a depth of 28 m in a homogeneous layer of London clay, the response is not affected by resonance in the soil. The wave pattern is similar to that caused by a point load in the soil. The propagation of waves into the soil is spherical but not symmetric, due to the eccentricity of the load on the tunnel invert with respect to the centre of the tunnel.
<figure><img src="image_10.png" /><figcaption>Fig. 9. Displacements of the Bakerloo line tunnel and the soil due to a harmonic excitation on the tunnel invert at 10 Hz for (a) $t = 0$ and (b) $t = 2T/8$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Displacements of the Bakerloo line tunnel and the soil due to a harmonic excitation on the tunnel invert at 40 Hz for (a) $t = 0$ and (b) $t = 2T/8$.</figcaption></figure>

### 3.6. Floating slab track

Next, the vibration isolation efficiency of a floating slab track with a resonance frequency of 11.68 Hz installed in the Bakerloo line tunnel is considered, using the Craig-Bampton substructuring technique. This example should be regarded as illustrative rather than practical, as the small diameter of the tunnel does not allow to install a floating slab track.

Figs. 11 and 12 show the displacements of the tunnel and the soil for a harmonic excitation at 10 and 40 Hz after installation of a floating slab in the tunnel. It should be noted that the grey scales are different from that in Figs. 9 and 10, as the displacements are generally lower. The floating slab track distributes vibrations along the tunnel that are radiated into the soil in directions parallel to the tunnel. Fig. 13 compares the displacements at the tunnel apex and in the free field above the tunnel for the cases without and with floating slab track; the isolation effect above 20 Hz can clearly be observed.

## 4. Numerical results for the RER B tunnel in Paris

### 4.1. Site characteristics

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 surface and a width of 11.9 m (Fig. 14). The
<figure><img src="image_12.png" /><figcaption>Fig. 11. Displacements of the Bakerloo line tunnel with floating slab track and the soil due to a harmonic excitation on the tunnel invert at 10 Hz for (a) $t = 0$ and (b) $t = 2T/8$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Displacements of the Bakerloo line tunnel with floating slab track and the soil due to a harmonic excitation on the tunnel invert at 40 Hz for (a) $t = 0$ and (b) $t = 2T/8$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Displacements (a) at the tunnel apex and (b) at the free field above the tunnel for the Bakerloo line tunnel without (solid line) and with (dashed line) floating slab track.</figcaption></figure>

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 = 14,000$ MPa, a Poisson's ratio $v^t = 0.15$, a density $\rho^t = 2400$ kg/m$^3$ and a hysteresis material damping ratio $\beta^t = 0.02$. The track is a classical ballasted track with UIC 60 rails supported every 0.60 m by grooved rubber pads on monobloc concrete sleepers.
<figure><img src="image_15.png" /><figcaption>Fig. 14. Cross-section of the metro tunnel on the line RER B of RATP at Cité Universitaire.</figcaption></figure>

The tunnel is embedded in a shallow layer of fill material (thickness $d = 1.6$ m), a layer of Beauchamp sand ($d = 3.2$ m) and a stiffer layer with marl and gravel ($d = 7.8$ m) on top of chalk. A SASW test revealed a shallow layer with $d = 1.4$ m and a shear wave velocity $C_s = 115$ m/s on top of a layer with $d = 2.8$ m and $C_s = 220$ m/s on a half-space with $C_s = 315$ m/s [38]. A Poisson's ratio $v^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. More recent down-hole measurements have confirmed the variation of the shear wave velocity with depth [39], and have demonstrated that the longitudinal wave velocity increases from $700$ m/s at a depth of $3$ m to $900$ m/s at a depth of $10$ m, resulting in a value of Poisson's ratio that decreases from $0.43$ to $0.33$ (and showing no evidence of the presence of ground water). A higher average density of $1902$ kg/m$^3$ was determined from soil samples, but has not been used in the calculations.

## 4.2. Kinematics of the tunnel

The cut-and-cover masonry tunnel is invariant in the $y$-direction. The periodicity introduced by the discrete support of the rails is neglected, since the track is not included in the model. The tunnel is invariant in the $y$-direction, so that the period $L$ can arbitrarily be chosen as $0.3$ m; the influence of the assumed period $L$ of the tunnel is presently investigated by a comparison with the results of a 2.5D calculation based on the semi-analytical insertion loss model [9–11]. The FE model of the tunnel consists of eight-node isoparametric brick elements. As no incompatible modes are included in the formulation of these elements, three elements are used over the height of the walls and ceilings, and two elements for the bottom slab, in order to represent the bending behaviour. The size of the finite elements is governed by the boundary element mesh along the tunnel-soil interface, so that minimum eight elements are used per wavelength at a maximum frequency of $80$ Hz; the shear wave velocity $C_s = 220$ m/s in the soil layer around the apex of the tunnel has been used to compute this wavelength.

Fig. 15 shows the first two in-plane and the first out-of-plane flexible modes $\Psi_m^0(\bar{x})$ of the reference cell. The response of the tunnel-soil system due to a harmonic load on the tunnel invert (in the centre of the left track in Fig. 14) will be computed in the frequency range upto $80$ Hz. Convergence analysis has shown that $30$ tunnel modes need to be included in the analysis [22].

## 4.3. Impedance of the tunnel

Fig. 16a shows the modulus of the sixth diagonal element $Z_t(6, 6)$ of the tunnel impedance matrix $Z_t(\kappa, \omega) = K_t(\kappa) - \omega^2 M_t(\kappa)$ as a function of the frequency and the slowness $p$. The element $Z_t(6, 6)$ corresponds to the
<figure><img src="image_16.png" /><figcaption>Fig. 15. The first (a,b) in-plane and (c) out-of-plane modes of the reference cell of the RER B tunnel. (a) Mode 5 at 9.90 Hz, (b) mode 6 at 13.96 Hz and (c) mode 9 at 41.95 Hz.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16. Modulus of the element Z<sub>t</sub>(6, 6) of the tunnel impedance matrix of the RER B tunnel (a) as a function of ω and p and (b) as a function of ω for p = 0 s/m (solid line), p = 0.4 × 10<sup>-3</sup> s/m (dashed line) and p = 0.7 × 10<sup>-3</sup> s/m (dotted line). Both figures are on a logarithmic scale.</figcaption></figure>

sixth mode of the reference cell (Fig. 15b) involving in-plane bending of the slab of the tunnel. Fig. 16a clearly illustrates that the sixth eigenfrequency of the tunnel reference cell corresponds to the cut-on frequency of a dispersive mode that propagates in the longitudinal direction of the tunnel. The phase velocity is infinite at the cut-on frequency and tends to a constant value of 1/C<sub>s</sub> = 0.628 × 10<sup>-3</sup> s/m for limiting high values of the frequency, where C<sub>s</sub> is the shear wave velocity in the masonry. Fig. 16b shows the variation of the element Z<sub>t</sub>(6, 6) of the tunnel impedance matrix with the frequency ω for three values of the slowness p. Similar figures can be made for other modes of the tunnel's reference cell, but are not shown here.

Fig. 17 shows, on a logarithmic scale, the modulus of the vertical displacement of the apex of the tunnel as a function of ω and p for the free tunnel that is not supported by the soil, including the
<figure><img src="image_18.png" /><figcaption>Fig. 17. Modulus of the vertical displacement of the apex of the RER B tunnel as a function of $\omega$ and $p$, accounting for 30 tunnel modes.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 18. (a) Real and (b) imaginary parts of the element $K_s(6,6)$ of the soil impedance matrix of the RER B tunnel as a function of $\omega$ and $p$.</figcaption></figure>

contribution of 30 modes of the tunnel reference cell. The response is governed by dispersive tunnel modes, of which the slowness tends to a value $1/C_s = 0.628 \times 10^{-3}$ s/m for the in-plane modes and $1/C_p = 0.403 \times 10^{-3}$ s/m for the out-of-plane modes, with $C_s$ and $C_p$ the shear and longitudinal wave velocities in the masonry.

## 4.4. Impedance of the soil

Fig. 18 shows the real and imaginary parts of the element $K_s(6,6)$ of the impedance matrix $K_s(\kappa, \omega)$ of the soil, corresponding to the sixth mode of the RER B tunnel. The real part increases for increasing values of $p$. High negative values of the imaginary part reflect dissipation of energy due to radiation damping or wave propagation in the soil, while the low imaginary part for high values of the slowness $p$ (or low values of the phase velocity along the tunnel axis) represents the absence of this radiation damping. As the RER B tunnel is shallow and embedded in a layered soil medium, the behaviour of the soil impedance as a function of $\omega$ and $p$ is more complicated than for the Bakerloo line tunnel, which was a deep tunnel embedded in a homogeneous soil layer. The soil impedance is now not only affected by the heterogeneous soil around the tunnel, but also by the proximity of the free surface.
<figure><img src="image_20.png" /><figcaption>Fig. 19. Transfer functions at the points (0, 0, 0) (bulleted line), (−2, 0, 0) (solid line), (−5, 0, 0) (dotted line), (−10, 0, 0) (dashed line) and (−20, 0, 0) (dashed-dotted line) along the free surface.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 20. Transfer functions at the points (−2.50, −0.15, −8.25) on the tunnel invert (bulleted line), (0, −0.15, −2.30) on the tunnel apex (dashed line) and (0, 0, 0) on the free surface (solid line).</figcaption></figure>

### 4.5. Response due to harmonic loading

Fig. 19 shows the transfer functions (vertical displacements) at some points along the free surface in the cross-section where the load is applied, while Fig. 20 shows the transfer functions at points on the tunnel invert, the tunnel apex and the free surface. The peak at 14 Hz is due to resonance of the soil above the tunnel; Fig. 20 confirms that the displacements of the tunnel apex and the soil above the tunnel are equal. Fig. 21 shows the displacements of the tunnel and the soil due to a harmonic loading on the tunnel invert at 14 Hz. The soil above the tunnel can be considered as a mass, which is moving in phase with the tunnel, that behaves as a spring. The peak at 14 Hz corresponds to the eigenfrequency of this equivalent single degree of freedom system (SDOF) system, while damping is due to the radiation of waves away from the tunnel.

Fig. 22 shows similar results at a frequency of 80 Hz. On the free surface above the tunnel, higher phase velocities (about 1920 m/s) are observed along the tunnel axis than in the direction perpendicular to the tunnel (about 650 m/s), due to the bending stiffness of the tunnel, resulting in an elliptical wave front. At a frequency of 80 Hz, the displacement of the ground surface and the tunnel apex are almost in phase, as this frequency
<figure><img src="image_22.png" /><figcaption>Fig. 21. Displacements of the RER B tunnel and the soil due to a harmonic excitation on the tunnel invert at 14 Hz for (a) $t = 0$ and (b) $t = 2T/8$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 22. Displacements of the RER B tunnel and the soil due to a harmonic excitation on the tunnel invert at 80 Hz for (a) $t = 0$ and (b) $t = 2T/8$.</figcaption></figure>

almost corresponds to the third vertical mode of the soil above the tunnel, as has been demonstrated by analysis of the dispersive characteristics of the waves in the soil above the tunnel [40].

## 5. Conclusion

A periodic coupled finite element-boundary element formulation is used to study the dynamic interaction between a tunnel and a layered soil due to a harmonic excitation on the tunnel invert. Two cases have been considered: a deep-bored tunnel on the Bakerloo line of London Underground and a shallow cut-and-cover masonry tunnel on the RER B line of RATP in Paris.

The Bakerloo line tunnel has a complex structure as the cast iron lining is stiffened with longitudinal and circumferential beams, resulting in a periodic structure. The periodic finite element formulation offers full modelling flexibility on a reference cell of the tunnel, which is an important advantage with respect to full three-dimensional or invariant tunnel models, which can be considered as the two extremes of modelling complexity. The soil impedance is formulated with a boundary element formulation on the soil-tunnel interface in the reference cell, allowing to model the soil as a horizontally layered half-space. This is still quite demanding from the computational point of view, however, depending on the size of the boundary element mesh and the evaluation of the Green-Floquet functions. As the Bakerloo line tunnel is an embedded tunnel, the interpretation of the soil's impedance benefits from the analytical study of the response of a homogeneous elastic medium with a cylindrical cavity. Analysis in the frequency-slowness domain demonstrates that, for the case under consideration, the soil dominates the response at low values of the slowness, while the tunnel dominates the response at high values of the slowness. The resulting response of the coupled soil-tunnel system is not dominated by dispersive modes in the tunnel and the soil. This affects the results in the frequency-spatial domain, showing spherical wave patterns in the soil, as would also result from a point load applied to the soil. These preliminary conclusions depend of course on the relative stiffness of the tunnel and the soil, as well as other factors such as the amount of material and radiation damping in the soil. The general
expectation is that, at low frequencies, the influence of in-plane deformations of the tunnel is low, whereas more complicated mode shapes are excited at higher frequencies, affecting the wave field radiated into the soil. However, this could not be clearly observed on the results for the Bakerloo line tunnel.

The case of the RER B tunnel in Paris is more complicated, as a shallow cut-and-cover tunnel with a non-circular cross-section is embedded in a heterogeneous layered soil. The only simplification here is that the tunnel is invariant in the longitudinal direction, allowing more freedom in the choice of a (shorter) period (which has not been studied in this paper and constitutes a subject of further research). It has been demonstrated that the response of this shallow tunnel is importantly influenced by the proximity of the free surface, giving rise to resonance phenomena, as well as to elliptical wave fronts radiated from the tunnel that originate from different phase velocities in the longitudinal and perpendicular direction.

It has also been demonstrated with an example how the Craig–Bampton substructuring technique allows to efficiently incorporate a track in the tunnel and to analyse the vibration isolation efficiency of a floating slab track in the tunnel, without the need to recompute the impedance of the soil. This will substantially facilitate future parametric studies on the efficiency of vibration isolation measures in the tunnel.

## 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 kindly acknowledged.

## References

[1] K.H. Chua, T. Balendra, K.W. Lo, Groundborne vibrations due to trains in tunnels, *Earthquake Engineering and Structural Dynamics* 21 (5) (1992) 445–460.

[2] K.H. Chua, K.W. Lo, T. Balendra, Building response due to sub-way train traffic, *Journal of Geotechnical Engineering, Proceedings of the ASCE* 121 (11) (1995) 747–754.

[3] W. Rücker, S. Said, *Erschütterungsübertragung zwischen U-Bahn-Tunneln und dicht benachbarten Gebäuden*, Forschungsbericht 199, Bundesanstalt für Materialforschung und -prüfung, Berlin, 1994.

[4] O. Von Estorff, A.A. Stamos, D.E. Beskos, H. Antes, Dynamic interaction effects in underground railway traffic systems, *Earthquake Engineering and Structural Dynamics* 8 (4) (1991) 167–175.

[5] C.J.C. Jones, D.J. Thompson, M. Petyt, Studies using a combined finite element and boundary element model for vibration propagation from railway tunnels, in: G. Guidati, H. Hunt, H. Heller, A. Heiss (Eds.), *Seventh International Congress on Sound and Vibration*, Garmisch-Partenkirchen, Germany, July 2000.

[6] A.A. Stamos, D.E. Beskos, Dynamic analysis of large 3-D underground structures by the BEM, *Earthquake Engineering and Structural Dynamics* (1995) 1–18.

[7] A.A. Stamos, D.E. Beskos, 3-D seismic response analysis of long lined tunnels in half-space, *Soil Dynamics and Earthquake Engineering* 15 (1996) 111–118.

[8] J.A. Forrest, Modelling of Ground Vibration from Underground Railways, Ph.D. Thesis, Department of Engineering, University of Cambridge, 1999.

[9] M.F.M. Hussein, H.E.M. Hunt, An insertion loss model for evaluating the performance of floating-slab track for underground railway tunnels, in: *10th International Congress on Sound and Vibration*, Stockholm, Sweden, July 2003 (CD-ROM).

[10] M.F.M. Hussein, H.E.M. Hunt, Dynamic effect of slab discontinuity on underground moving trains, in: *11th International Congress on Sound and Vibration*, St. Petersburg, Russia, July 2004, pp. 3045–3054 (CD-ROM).

[11] M.F.M. Hussein, Vibration from Underground Railways. Ph.D. Thesis, Department of Engineering, University of Cambridge, 2005.

[12] X. Sheng, C.J.C. Jones, M. Petyt, Ground vibration generated by a harmonic load acting on a railway track, *Journal of Sound and Vibration* 225 (1) (1999) 3–28.

[13] X. Sheng, C.J.C. Jones, M. Petyt, Ground vibration generated by a load moving along a railway track, *Journal of Sound and Vibration* 228 (1) (1999) 129–156.

[14] C. Madshus, A.M. Kaynia, High-speed railway lines on soft ground: dynamic behaviour at critical train speed, *Journal of Sound and Vibration* 231 (3) (2000) 689–701.

[15] G. Lombaert, G. Degrande, S. François, J. Kogut, L. Pyl, Validation of a numerical model for the prediction of railway induced vibrations in the free field, in: *11th International Conference on Soil Dynamics and Earthquake Engineering and Third International Conference on Earthquake Geotechnical Engineering*, Berkeley, CA, USA, January 2004.

[16] G. Lombaert, G. Degrande, J. Kogut, S. François, W. Haegeman, L. Karl, An experimental validation of a numerical model for railway induced vibrations at different levels, in: D. Thompson, C. Jones (Eds.), *Eighth International Workshop on Railway Noise, Vol. 1, Buxton, UK, September 2004*, pp. 97–106.
[17] G. Lombaert, G. Degrande, B. Vanhauwere, B. Vandeborgh, A numerical study of vibration isolation for rail traffic, in: *ISMA2004 International Conference on Noise and Vibration Engineering*, Leuven, Belgium, September 2004.

[18] D. Clouteau, D. Aubry, M.L. Elhabre, E. Savin, Periodic and stochastic BEM for large structures embedded in an elastic half-space, in: *Mathematical Aspects of Boundary Element Methods*, CRC Press, London, 1999, pp. 91–102.

[19] 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.

[20] D. Clouteau, M.L. Elhabre, D. Aubry, Periodic BEM and FEM-BEM coupling: application to seismic behaviour of very long structures, *Computational Mechanics* 25 (2000) 567–577.

[21] (http://www convurt com), 2003.

[22] D. Clouteau, M. Arnst, T.M. Al-Hussaini, G. Degrande, Free field vibrations due to dynamic loading on a tunnel embedded in a stratified medium, *Journal of Sound and Vibration* 283 (1–2) (2005) 173–199.

[23] D. Aubry, D. Clouteau, G. Bonnet, Modelling of wave propagation due to fixed or mobile dynamic sources, in: N. Chouw, G. Schmid (Eds.), *Workshop Wave '94, Wave propagation and Reduction of Vibrations*, Ruhr Universität Bochum, Germany, December 1994, pp. 109–121.

[24] P. Chatterjee, G. Degrande, Free field vibration measurements due to a test train at the site of Cité Universitaire in Paris, Report BWM-2002-09, Department of Civil Engineering, K.U. Leuven, November 2002 (CONVURT EC-Growth Project G3RD-CT-2000-00381).

[25] S. Jacobs, P. Chatterjee, G. Degrande, Vibrations due to service trains in the 'Maison du Mexique' on the RER B line at Cité Universitaire in Paris, Report BWM-2003-02, Department of Civil Engineering, K.U. Leuven, February 2003 (CONVURT EC-Growth Project G3RD-CT-2000-00381).

[26] P. Chatterjee, G. Degrande, S. Jacobs, J. Charlier, P. Bouvet, D. Brassenx, Experimental results of free field and structural vibrations due to underground railway traffic, in: *10th International Congress on Sound and Vibration*, Stockholm, Sweden, July 2003 (CD-ROM).

[27] G. Degrande, P. Chatterjee, S. Jacobs, J. Charlier, P. Bouvet, D. Brassenx, CONVURT project—experimental results of free field and structural vibrations due to underground railway traffic, in: *World Congress on Railway Research*, Edinburgh, UK, September 2003.

[28] P. Chatterjee, G. Degrande, S. Jacobs, Free field and building vibrations due to the passage of test trains at the site of Regent's Park in London, Report BWM-2003-20, Department of Civil Engineering, K.U. Leuven, December 2003 (CONVURT EC-Growth Project G3RD-CT-2000-00381).

[29] G. Degrande, P. Chatterjee, W. Van de Velde, P. Hölscher, V. Hopman, A. Wang, N. Dadkah, R. Klein, Vibrations due to a test train at variable speeds in a deep bored tunnel embedded in London clay, in: *11th International Congress on Sound and Vibration*, St-Petersburg, Russia, July 2004, pp. 3055–3062 (CD-ROM).

[30] R.J. Craig, M. Bampton, Coupling of substructures for dynamic analyses, *AIAA Journal* 6 (7) (1968) 1313–1319.

[31] D. Clouteau, G. Degrande, G. Lombaert, Numerical modelling of traffic induced vibrations, *Meccanica* 36 (4) (2001) 401–420.

[32] X. Sheng, C.J.C. Jones, D.J. Thompson, Responses of infinite periodic structures to moving or stationary harmonic loads, *Journal of Sound and Vibration* 282 (1–2) (2005) 125–129.

[33] A.K. Chopra, *Dynamics of Structures: Theory and Applications to Earthquake Engineering*, Prentice-Hall, Englewood Cliffs, NJ, 1995.

[34] P. Hölscher, V. Hopman, Test site Regent's Park London, Soil description, Report 381540-104, Version 2, GeoDelft, December 2003 (CONVURT EC-Growth Project G3RD-CT-2000-00381).

[35] L. Pyl, G. Degrande, Determination of the dynamic soil characteristics with the SASW method at Regent's Park in London, Report BWM-2003-17, Department of Civil Engineering, K.U. Leuven, December 2003 (CONVURT EC-Growth Project G3RD-CT-2000-00381).

[36] K.J. Bathe, *Finite Element Procedures*, Prentice-Hall, Englewood Cliffs, NJ, 1996.

[37] A. Boström, A. Burden, Propagation of elastic surface waves along a cylindrical cavity and their excitation by a point force, *Journal of the Acoustical Society of America* 72 (3) (1982) 998–1004.

[38] L. Pyl, 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, K.U. Leuven, October 2002 (CONVURT EC-Growth Project G3RD-CT-2000-00381).

[39] P. Hölscher, V. Hopman, Test site Cité Universitaire Paris, Soil description, Report 381540-0132, Draft version, GeoDelft, June 2004 (CONVURT EC-Growth Project G3RD-CT-2000-00381).

[40] M. Arnst, Three-dimensional modelling of free field and structural vibration due to harmonic and transient loading in a tunnel, Master's Thesis, Department of Civil Engineering, K.U. Leuven, Erasmus, Ecole Centrale de Paris, France, 2003.