Modelling ground vibration from railways using wavenumber finite- and boundary-element methods

X Sheng, C.J.C Jones, D.J Thompson


A mathematical model is presented for ground vibration induced by trains, which uses wavenumber finite- and boundary-element methods. The track, tunnel and ground are assumed homogeneous and infinitely long in the track direction (x-direction). The models are formulated in terms of the wavenumber in the x-direction and discretization in the yz-plane. The effect of load motion in the x-direction is included. Compared with a conventional, three-dimensional finite- or boundary-element model, this is computationally faster and requires far less memory, even though calculations must be performed for a series of discrete wavenumbers. Thus it becomes practicable to carry out investigative study of train-induced ground vibration. The boundary-element implementation uses a variable transformation to solve the well-known problem of strongly singular integrals in the formulation. A ‘boundary truncation element’ greatly improves accuracy where the infinite surface of the ground is truncated in the boundary-element discretization. Predictions of vibration response on the ground surface due to a unit force applied at the track are performed for two railway tunnels. The results show a substantial difference in the environmental vibration that could be expected from the alternative designs. The effect of a moving load is demonstrated in a surface vibration example in which vibration propagates from an embankment into layered ground.


1. Introduction

An important issue with regard to the environmental impact of rail traffic is the problem of ground vibration produced by trains running in tunnels or on the ground surface. The dominant frequencies cover a range of about 2 to 200 Hz (Grootenhuis 1977). Components of low frequencies in this range cause whole-body vibration, while those of higher frequencies give rise to structure-borne noise. Experimental evaluations of vibration-reducing designs are notoriously difficult to obtain since ground and construction parameters are difficult to control. To understand the mechanisms involved in vibration generation and propagation, and to develop well-founded vibration mitigation solutions, a realistic mathematical model is needed.

For trains running on the surface of a layered ground, models have been developed for ground vibration generated by moving axle loads (Krylov 1995, 1996; Degrande & Lombaert 2001), moving dynamic loads (Sheng et al. 1999) and by both the moving axle loads and wheel and rail irregularities (Sheng et al. 2003a). These models employ a semi-analytical approach in the wavenumber-frequency domain.

Semi-analytical models have also been derived for trains running in tunnels. Forrest (1999) presents a model based on cylinder theory, in which the ground and the tunnel were modelled as a viscoelastic whole space including a circular thin shell of infinite length. Such an analytical model is computationally very efficient but its usefulness is limited if the tunnel has a non-circular cross-section, if the tunnel is not well below the ground surface or if the ground is layered. Another model (Metrikine & Vrouwenvelder 2000) is composed of two two-dimensional viscoelastic layers, representing the soil above and below the tunnel, and two identical beams representing the tunnel. The beams are connected by continuously distributed springs. This model is able to account for the movement of excitations and wave propagation in the tunnel direction, but not for the wave propagation in the lateral direction.

For ground vibration from tunnels, numerical approaches are normally employed to allow arbitrary geometry. Most commonly, either the finite-element method (FEM) (Gardien & Stuit 2003) or the boundary-element method (BEM) (Andersen & Jones 2002) is used.

In FEM, built structures and part of the ground are discretized, introducing an artificial boundary. Proper boundary conditions, often termed ‘transmitting boundary conditions’, are required in order to ensure that no wave of significant amplitude is reflected. In general, however, exact transmitting boundary conditions cannot be found. Approximate boundary conditions may be used but they require that the artificial boundary be placed in the far field, which increases the size and complexity of the analysis. For a three-dimensional analysis of a large structure, the number of elements involved in the analysis rapidly becomes too large as frequency increases.

The BEM for the analysis of solid elastic bodies has been well set out in recent years (Dominguez 1993; Gaul et al. 2003) and is very well suited to the analysis of infinite media, since the radiation of waves towards infinity is automatically included. However, BEM is not optimal for thin structures such as tunnels since both faces of a structure have to be discretized and BEM has numerical problems for thin domains. Therefore, for the dynamics of soil–structure interaction, a combination of FEM and BEM is a natural choice in which the finite structure is modelled using FEM and the ground is modelled using BEM. A three-dimensional finite-element and boundary-element (FE–BE) model has been implemented by Andersen & Jones (2002) for the study of train-induced ground vibration. However, this model is still very complicated and expensive to use, and is therefore impractical for carrying out investigative studies for the whole frequency range of interest.

In many practical situations, the geometry of the ground-tunnel structure can be assumed to be invariant with respect to translation along the tunnel axis. This feature is also present for track structures (rails, ballast and embankment) and for many vibration mitigation devices, such as ditches and wave impeding blocks (Takemiya 2002). For such structures, an approach may be applied in which the problem is transformed into a sequence of two-dimensional (yz-plane) models, depending on the wavenumber in the tunnel direction (x-direction). For each of a series of discrete wavenumbers in the x-direction, the finite cross-section of the built structure is modelled using the FEM and the wave propagation in the surrounding soil is modelled using the BEM. The global FE–BE equations are coupled and then solved, giving the component of the response at that wavenumber. The responses at different x-positions are constructed from these components using an inverse Fourier transform scheme. These methods may be termed the wavenumber FE–BE methods or two-and-a-half dimensional FE–BE methods. As the discretization is only made over the cross-section of the ground-tunnel structure, the total number of degrees of freedom is greatly reduced compared with the corresponding three-dimensional model.

It should be noted that similar terms such as the discrete wavenumber method, the discrete wavenumber boundary-element method and the discrete wavenumber boundary integral equation method have also been used in the literature (Kim & Papageorgiou 1993; Haartsen et al. 1994; Bouchon 2003) to refer to the conventional BEM in which the half-space (homogeneous or layered) Green's functions are evaluated using the discrete wavenumber method proposed by Bouchon (2003).

An alternative scheme to the one presented in this paper was developed by the present authors (Sheng et al. 2003b) based on the ‘fictitious forces method’ of Luco & deBarros (1993), in which the Green's function for a layered half-space was used and the forces around the tunnel–soil interface were expressed as a series of harmonic terms rather than discretized in the finite-element sense. Because of the time taken to compute the layered half-space Green's function, the present method has been found to be more efficient computationally in example cases. It is also a simpler method within which to generalize the geometry.

The wavenumber method for FE was first proposed by Gavric (1994, 1995) for computing the dispersion curves of thin-walled waveguides and free rails. Recently, Shah et al. (2001) employed this method to analyse the wave propagation property of thin-walled structural members using in-plane and out-of-plane plate elements.

The formulation of the direct wavenumber BEM and its applications are mostly found in papers concerning seismic responses of scatterers with a two-dimensional topology subject to incident seismic plane waves. Stamos & Beskos (1996) consider only harmonic plane waves as both incident and scattered waves. The Fourier-transformed (with respect to x only) whole-space Green's functions are required owing to unit stationary harmonic forces and at a single wavenumber. That single wavenumber is determined by the frequency, speed and propagation direction of the incident wave. Owing to the translational invariance of the waves with respect to the x-coordinate, Papageorgiou & Pei (1998) consider each harmonic component of the scattered wave as being produced by a constant load moving uniformly in the x-direction at the incident wave speed projected in the same direction. Therefore, the Fourier-transformed whole-space Green's functions corresponding to a moving constant load are utilised in the boundary integral equation.

It is clear that a systematic derivation of the wavenumber BEM is needed for its application to train-induced ground vibration. As mentioned above, the Green's functions play a key role in the BEM. The Green's functions are available for a homogeneous full-space excited by stationary harmonic sources (Eason et al. 1955) and for a layered half-space subject to both stationary and moving harmonic sources (de Barros & Luco 1994; Sheng 2001; Ronald & Bojan 2002). The so-called two-and-a-half-dimensional Green's functions for a homogeneous whole-space have also been derived for stationary harmonic sources (Stamos & Beskos 1996; Tadeu & Kausel 2000) and recently for moving harmonic sources (Sheng et al. 2002).

This paper describes a model of ground vibration from tunnels using the wavenumber FE and (direct) BE methods. In §2, the wavenumber FE formulation is outlined briefly. A systematic formulation of the wavenumber BEM is presented in §3. Coupling between FE domains and BE domains is dealt with in §4. In §5, two tunnel designs are compared in terms of vibration on the ground surface using this model. Section 6 outlines an example of the effect of a moving load.

2. Wavenumber finite-element method

(a) Differential equation of motion of an element

Suppose that an elastic body is infinitely long in the x-direction and that its cross-sections normal to the x-axis are invariant. The x cross-section, A, is discretized into a number of elements. The same discretization is also made on the x+dx cross-section. An element on the x cross-section and its counterpart on the x+dx cross-section define an element prism. The displacements of the n nodes on the element are denoted by a 3n vector:Embedded Image(2.1)in which the first subscript denotes the component of displacement in the x-, y- and z-directions and the second denotes the node number.

The corresponding nodal force vector is denoted by {F(x,t)}. A shape function matrix of order 3×3n is defined and denoted by [Φ(y, z)], so that the displacement vector of the element at any point (x, y, z) may be approximated asEmbedded Image(2.2)

The kinetic energy, T, of the element prism (of length dx) is expressed in terms of the nodal displacements:Embedded Image(2.3)whereEmbedded Image(2.4)(ρ is the material density). If the stress–strain relationship in terms of E, the Young's modulus and ν, the Poisson's ratio, is given by Petyt (1990)Embedded Image(2.5)then the potential energy, U, is given byEmbedded Image(2.6)

whereEmbedded Image

The vector of generalized forces is derived by considering the virtual work done by the stresses on the x and x+dx cross-sections of the element prism.Embedded Image(2.7)whereEmbedded ImageThis leads to the expression of the generalized forces asEmbedded Image(2.8)

By inserting the expressions for T, U and Q into Lagrange's equation,Embedded Image(2.9)the differential equation of motion of the element is obtained:Embedded Image(2.10)in which the terms have been collected together so thatEmbedded Image(2.11)[M], [K]0 and [K]2 are 3n×3n symmetric matrices, [M] and [K]2 are positive definite and [K]0 is non-negative. It can also be shown that [K]1 is an antisymmetric matrix.

(b) Global finite-element equation

The conventional ‘assembly’ of the element matrices in equation (2.10) can be applied to obtain the corresponding matrices of the assembled FEM and thus, the global differential equation of motion. This is still represented by equation (2.10). The sub-matrices corresponding to y- and z-coordinates in the matrices [K]0 and [M] are those of the plane-strain problem.

Now it is assumed that the nodal forces apply at points moving at constant speed c in the x-direction (i.e. it is implied that the nodes of the model represent such moving points) and that, in the moving frame of reference, these forces are time-harmonic with radian frequency Ω. The Fourier transform of any quantity, f(x,t), at a node due to such loads is denoted by Embedded Image, that is,Embedded Image(2.12)where β is the wavenumber in the x-direction. It is shown in Dieterman & Metrikine (1997) and also Sheng et al. (1999) that Embedded Image is harmonic with a modified frequency ω=Ωβc. In other words, it can be written thatEmbedded Image(2.13)where Embedded Image is the amplitude of Embedded Image. In what follows, Embedded Image is also termed the Fourier transform of f(x,t), though the factor Embedded Image has been dropped. Embedded Image is used to represent the correspondingly transformed displacements q(x, t). Substitution of equations (2.12) and (2.13) into equation (2.10), and adopting the bold-italic letter notation for global vectors and matrices in the frequency–wavenumber domain (i.e. the domain of assembly and solution of system matrices) yieldsEmbedded Image(2.14)

The transformed displacement vector Embedded Image can be evaluated from equation (2.14) for each β and then the actual displacements may be obtained using an inverse Fourier transform.

3. Wavenumber boundary-element method

(a) The two-and-a-half-dimensional reciprocal theorem in elasto-dynamics

As in §2, the x cross-section of the infinitely long elastic body is denoted by A. Now, the boundary of A in its own plane is denoted by Γ.

For this elastic body, two elasto-dynamic states are defined. The first one is described by displacements uk(x,y,z,t), body forces ρbk(x,y,z,t) and boundary tractions pk(x,y,z,t), where k=1, 2, 3 corresponding to x-, y- and z-directions. The second state is described by Embedded Image, Embedded Image and Embedded Image. A reciprocal relation between these two states exists and is stated asEmbedded Image(3.1)(Archenbach 1973; Dominguez 1993), where the summation convention over the repeated index applies, * denotes the Riemann convolution which applies only with time t,Embedded Image(3.2)and the initial displacements and velocities areEmbedded Image(3.3)

Since the elastic body is infinitely long in the x-direction, the two elasto-dynamic states can be expressed in terms of the Fourier integral, for instance,Embedded Image(3.4)

Using this in the reciprocal relation, equation (3.1), givesEmbedded Image(3.5)

The triple integrals in equation (3.5) can be simplified by making use of the fact that Embedded Image and Embedded Image are independent of x, for example,Embedded Image(3.6)

SinceEmbedded Imagewhere δ(.) is the Dirac-δ and Embedded Image, Embedded Image, equation (3.6) becomesEmbedded Image(3.7)By inserting the first form of this into the left-hand side of equation (3.5) and the second form into the right-hand side, it is shown that equation (3.5) holds if, for every value of β,Embedded Image(3.8)

Now it is assumed that the first elasto-dynamic state is due to harmonic loads of radian frequency Ω moving in the positive x-direction at speed c, while the second state is due to harmonic loads of the same frequency but moving in the negative x-direction at speed c. Thus, as in equation (2.13), it may be writtenEmbedded Image(3.9)so that Embedded Image and Embedded Image. Inserting these expressions into equation (3.8) and performing the convolution defined by equation (3.2), yieldsEmbedded Image(3.10)This represents the reciprocal relation between the two states in the wavenumber domain. When β is set to zero, this recovers the reciprocal theorem in elasto-dynamics for the plane-strain problem.

(b) Integral representation

Now the displacements and tractions in the second elasto-dynamic state are chosen to be those of a whole-space due to a unit harmonic load of frequency Ω. This load moves at speed c in the negative x-direction along a straight line which passes through (y0,z0). Thus, the body force for the second state is given byEmbedded Image(3.11)where δkl is Kronecker's delta. The Fourier transform of equation (3.11) is given byEmbedded Image(3.12)

Insertion of equation (3.12) into equation (3.10) givesEmbedded Image(3.13)This expresses the displacements at point (y0, z0), in terms of the displacements and tractions on the boundary and the externally applied body forces. Embedded Image and Embedded Image are used to denote the Fourier-transformed (displacement and traction) Green's functions (Sheng et al. 2002) due to the moving load defined in equation (3.11); they are given in appendix A. The first subscript indicates the response direction while the second denotes the source direction.

Now a boundary integral equation is established by taking point (y0, z0) on to the boundary Γ so that equation (3.13) becomesEmbedded Image(3.14)where Embedded Image is the displacement vector at (y0, z0), Embedded Image is the surface traction vector, and the displacement and traction Green's function matrices areEmbedded Image(3.15)

(c) Discretization using quadratic shape functions

The boundary integral equation (3.14) is discretized using the three-noded quadratic boundary-element shown in figure 1. The shape function matrix is given by:Embedded Image(3.16)where [I] is the 3×3 unit matrix (if [Φ(ξ)] is used to approximate the boundary geometry, then [I] is the 2×2 unit matrix) andEmbedded Image(3.17)are shape functions (Dominguez 1993). The displacements, tractions and coordinates along a boundary-element may be approximated by those at the nodes of the element via the shape functions. By doing so, equation (3.14) becomesEmbedded Image(3.18)where NE is the number of elements, Embedded Image is the displacement vector of node i and Embedded Image and Embedded Image are two 9×1 vectors consisting of the nodal displacements and tractions, respectively, for the three nodes on element j. Embedded Image and Embedded Image are the displacement and traction Green's function matrices with the source position corresponding to node i (called the collocation point) and the observer position given byEmbedded Image(3.19)where {y}j is a 6×1 vector containing the coordinates of the element nodes. Finally,Embedded Image(3.20)represents the contribution to node i from the body forces.

Figure 1

Three-noded, quadratic shape function boundary-element.

For a system of elements with a total of N nodes, equation (3.18) is assembled for each collocation node i to form a global boundary-element equationEmbedded Image(3.21)where [H] and [G] are 3N×3N matrices.

(d) Evaluation of singular integral terms

In constructing the global equation (3.21), the integrals are evaluated for the Green's functions on each element making up the total integration around the boundary. For an element excluding the collocation point, the integration along this element can be carried out using a standard numerical quadrature. In the present study, a 10-point Gauss–Legendre quadrature rule is used, the high order reflecting the high polynomial order of the integrands.

When the element of integration contains the collocation node, the integrands in both [H] and [G] become singular. The singular terms for constructing [G] may be integrated after a special treatment that takes into account the presence in the integrands of a weak singularity which is of logarithmic order (Singh & Tanaka 2001; Johnston & Johnston 2002). The singular terms for constructing [H], however, are of higher order and are not integrable. The total integration around the boundary for these terms has a finite result because the integrands are either a δ-like function, or have singularities of opposite signs either side of the singular point and they cancel. The integral cannot be evaluated using any quadrature scheme on an element-by-element basis and an alternative method must be devised. For a stationary harmonic load, the difficulty may be overcome by employing the ‘rigid body motion technique’ combined with either a fully or a locally enclosing elements technique (Dominguez 1993; Andersen & Jones 2002). In this paper, use is made of a new and more efficient technique, which is described below.

The irregularities of tractions Embedded Image at and near the collocation point (i.e. as r→0) where the boundary is smooth are characterized asEmbedded Image(3.22)where p2 is defined in appendix A, a and b are regular at r=0, δij is the Kronecker delta and δ(.) is the Dirac delta-function. The traction expressions derived from the applied displacements via Hooke's law do not include the delta-function part. This part only arises in the case of an applied point force at the collocation point. In what follows, the traction expressions derived from the displacements are termed the base tractions. Two cases are dealt with separately below.

  1. When the collocation point is the middle node (node 2) of element Γ1, the 3×3 sub-matrices associating nodes p (1, 2 and 3) in matrix [H] are evaluated on a local coordinate system and, according to equation (3.18), are given byEmbedded Image(3.23)where Embedded Image is the Jacobian, that is, Embedded Image. Note that the unit load at the collocation point is treated as a traction rather than a body force, therefore the Embedded Image term in equation (3.18) must be dropped. As ξ→0, ϕ1(ξ)∼−0.5ξ, ϕ3(ξ)∼−0.5ξ (equation (3.17)), the contribution to the integrals in equation (3.23) for p=1 or 3 from the delta-function part in equation (3.22) becomes zero. Thus it is only possible to use the base tractions in equation (3.23); the resulting integrands are regular at ξ=0 and can be evaluated using a conventional quadrature rule. Since ϕ2(ξ)=ϕ2(−ξ), for p=2 equation (3.23) may be written asEmbedded Image(3.24)

    Now substitute Embedded Image into equation (3.24). For the delta-function part, the integral equals 0.5δij. For the base tractions, the integrand is irregular at ξ=0, at most of logarithmic order, due to the cancellation of the irregularities of order 1/r appearing in the two terms in the bracket. Thus, by employing the same transformation as used for the evaluation of the irregular integrals in matrix [G], the integral in equation (3.24) can be evaluated numerically without difficulty.

  2. When the collocation point is an end node of an element, say node 3 which joins element 1 and element 2, the elements in [H] associated with this node are given byEmbedded Image(3.25)since ϕ1(ξ)=ϕ3(−ξ). Again, for the delta-function part of the tractions, this gives 0.5δij, and for the base tractions, the integrand in this equation is irregular at ξ=1, at most of the logarithmic order. Thus, it too can be evaluated numerically.

4. Coupling between subdomains, reduction of the truncation effect and inclusion of damping

(a) Coupling between subdomains

Up to this point, the BEM formulation has been derived for a domain consisting of a single homogeneous material. A domain such as a layered ground, which may also incorporate built structures such as a tunnel and a track, can be divided into a number of BE and FE subdomains. Each of the BE subdomains is homogeneous. For each subdomain, a BE equation or a FE equation can be constructed, and coupling of these equations gives the global equation for the whole system.

A single BE equation may readily be constructed for all the BE subdomains. From equation (3.21) pre-multiplied by [G]−1, this equation may be written in the formEmbedded Image(4.1)

Note that equilibrium at the interfaces between assembled BE domains causes corresponding elements of Embedded Image to vanish. If, in addition to the BE subdomains, an FE subdomain is present, then the FE equation must be coupled with equation (4.1) to give the global equation for the whole model. The displacements, tractions and nodal forces of the nodes at the FE–BE interface are denoted by Embedded Image, Embedded Image and Embedded Image. Those of the remaining nodes are denoted by Embedded Image and Embedded Image for the BE subdomains, and Embedded Image, Embedded Image for the FE subdomain. Thus equation (4.1) can be split intoEmbedded Image(4.2)and the FE equation from equation (2.14) intoEmbedded Image(4.3)

A transformation matrix, [T], may be constructed to convert the tractions Embedded Image of the boundary-element formulation at the FE–BE interface into the equivalent nodal forces to enable assembly with the FE matrices, that is,Embedded Image(4.4)

The dimension of the matrix [T] depends on the number of nodes at the FE–BE domain interface. Inserting equation (4.4) into equation (4.2) givesEmbedded Image(4.5)The tractions on the FE–BE interface have thus been converted into equivalent nodal forces. This means the possible discontinuity of tractions on the FE–BE interface does not have to be considered. Equations (4.3) and (4.5) yieldEmbedded Image(4.6)This is the global equation for the whole system.

In many situations, the track-tunnel-ground system, the loads applied and the response may be considered to be symmetric about the xz-plane. Reduction of the dimension (by about half) of the global equation (4.6) is achieved by making use of the symmetry, so that the computational efficiency can be greatly improved.

(b) Reduction of truncation effect

In the BE model, since the whole-space Green's functions are used, the ground surface and the horizontal interfaces of a layered ground are represented by boundary-elements and therefore must be truncated. The truncation may give significant errors to the calculated results near the edges of the model for a layered ground (Andersen & Jones 2002). A special element has therefore been developed in order to reduce the error due to the truncation.

Assume, as shown in figure 2, the nodes on the edge element are numbered as 1, 2 and 3. Beyond the edge element, extra elements may be implemented on the boundary (ground surface or interface) with nodes having numbers 0, −1, −2, …, to construct an extended boundary-element model. The boundary-element equations of the extended model for collocation points being nodes i are given byEmbedded Image(4.7)where i=1, 2, …, N and N is the number of the nodes on the original boundary-element model. Dropping the terms associated with the extra nodes, equation (4.7) becomes the boundary-element equation of the original model. Now a wave field is assumed for the extra elements (called the far field),Embedded Image(4.8)where yj>y1 and j=0, −1, −2, …. If the shear wave speed of the material is denoted by c2, then Embedded Image is the complex wavenumber in the y-direction, where ω=Ωβc as before. With the substitution of this and equation (4.8) into equation (4.7), modified BE equations for the original model are produced:Embedded Image(4.9)

Figure 2

Nodes on the edge of the BE.

The extent to which equation (4.9) is more accurate than the original BE equation depends on how well the far field wave is modelled by the assumed expressions in equation (4.8). In actual calculation, only the equations associated with the nodes on the edge element need modification, that is, the last element of the truncated boundary is replaced by a ‘truncation element’. Experience of a number of example calculations shows that by this treatment, the accuracy of the calculated results near the edges of a BE model can be greatly improved.

(c) Inclusion of damping

The effect of damping can be included in both the FEM and in the BEM as a hysteretic material model, that is, the Young's modulus is specified as Embedded Image where η is the material loss factor. Although other models of damping may be proposed, the imprecise character of the material parameters for soils (which is the application here) does not support the adoption of such a model; most researchers resort to the use of a loss factor as a practical descriptor valid for the frequency range of interest for environmental vibration from trains.

It is well known that this leads to an error in the form of a non-causal solution; a response of the structure to the load is found even before the load is applied (Crandall 1970). This presents a difficulty when a moving load is modelled. The magnitude of the non-causal response has been studied by Sheng (2001) for impulse responses of a one degree of freedom system model, a model of a railway track on a Winkler foundation, a whole-space elastic medium and an elastic half-space, in terms of the ratio of the power of the response before or after the impulse (or arrival of the impulse at a distance). It was concluded that for loss factors below about 0.2, the effect is sufficiently small to be neglected for practical model cases.

5. Results

(a) Validation against exact models

In order to check the accuracy of the method of combined FEM and BEM, two analyses have been performed for special cases that can be compared with exact calculations. This also provides a test of the element size necessary for sufficient accuracy at the highest frequency of interest, in this case 200 Hz.

In the first case, the response of a homogeneous half-space, subject to a stationary vertical point load at a depth of 18 m, is calculated. The half-space has a Young's modulus of 1.77×109 N m−2, a Poisson's ratio of 0.4 and a density of 1700 kg m−3, corresponding to a P-wave speed of 1500 m s−1 and S-wave speed of 610 m s−1. Material damping is included in the ground by use of a damping loss factor of 0.15.

In the BEM, the vertical symmetry about the load point is used and the surface of the ground is modelled to a distance of 50 m to the side of the load position using 50 elements. At this point an edge element is used to terminate the BEM of the surface. The element length used for all the elements corresponds to three elements per Rayleigh wavelength at 200 Hz. The response has been evaluated for 512 wavenumbers with a spacing of β of 0.005×2π rad m−1. The reverse FFT therefore yields the responses at 512 positions in x from −100 to +100 m.

The exact result for this case has been evaluated using the wavenumber-domain three-dimensional flexibility matrix method described by Sheng (2001). This model is ‘semi-analytical’; it too uses a reverse FFT to retrieve responses as a function of Cartesian coordinates from the analytical wavenumber-domain solution. The instantaneous vertical displacements on the ground surface (in phase with the load) along the x-axis and y-axis at 200 Hz are plotted in figure 3. The BE result is shown with and without the use of the edge termination element. This case shows the importance of the edge termination element for the results in the discretized yz-plane. The results in the xz-plane are not affected by the BE surface truncation. With the edge termination element, it can be seen that the BEM model has good accuracy, even at this frequency in both directions and even at the edge of the model.

Figure 3

Vertical displacement on the ground surface (a) along the x-axis, (b) along the y-axis. ---, using the semi-analytical method; ---, using the BEM with the edge termination element; ⋯, without the edge termination element.

In the second test case, the FEM part of the model is also tested. A homogeneous whole space is modelled with an infinitely long circular lined ‘tunnel’ (a cylindrical tube). In this case, the tube has inner and outer radii of 6.65 and 7.25 m and is modelled using 30 eight-noded (distorted) quadrilateral finite-elements, each approximately 1.45×0.6 m. The solution has been carried out for the same ground properties as in the first case. The tube has been ascribed a Young's modulus of 40×109 N m−2, Poisson's ratio of 0.15 and density of 2400 kg m−3, representing concrete.

Again, for comparison with the FEM–BEM model, a model that is analytical in the wavenumber domain is used and the response is recovered by reverse FFT (Sheng et al. 2002). The vertical displacements along the bottom of the tube resulting from a unit vertical point load of 200 Hz applied at the bottom of the tube at x=0 are shown in figure 4. This indicates good accuracy has been achieved by the FE–BE models.

Figure 4

Vertical displacement along the tunnel invert. —, using the semi-analytical method; - - -, using the FE–BE method (curves coincide except near x=0 m).

(b) Comparison of the vibration response from two tunnels

In order to demonstrate a more practical use of the modelling method, a comparative analysis is presented of two alternative tunnel designs. A single double-track tunnel and a pair of single-track tunnels have been considered as alternative designs for a double-track line carrying high speed vehicles. The cross-sections of the tunnel designs are shown in figure 5.

Figure 5

(a) The single bore tunnel design. (b) the twin bore tunnel design.

The inner and outer radii of the single bore tunnel ring are 6.65 and 7.25 m and those of the twin bores are 4.3 and 4.75 m. The axes of the two tunnels of the twin bore design are 21.5 m apart. Connecting passages have been ignored in the analysis of the twin bore case. The depth of the railhead below the ground surface is taken to be 28.52 m for both tunnels. Having the same depth at railhead level is the comparative situation that arises in practice, since the vertical alignment of the track is fixed. This means that the crowns of the twin bore tunnels are 21 m from the ground surface and the crown of the single bore tunnel is only 18 m from the surface. To compare the tunnels' vibration at the surface, the ground is modelled only as a simple half-space of soil with a shear wave speed of 700 m s−1, a dilatational wave speed of 1710 m s−1 and a mass density of 2000 kg m3. A damping loss factor of 0.05 is used to represent the material losses in the soil. The concrete in each of the structures has the same properties as those of the ‘tube’ in §5a.

For the twin bore design, 80 BEs are used for the ground surface from the vertical centre-line of the two tunnels to a distance of 80 m, 30 FEs for the tunnel ring and six FEs for the tunnel invert. A similar element size is used for the model of the single bore tunnel. Note that, since the tunnel and ground structures are symmetrical about the xz-plane, only half the structure has to be discretized. Despite using the symmetry in the geometry of the tunnel, the loading is not symmetric, that is, only the track in the discretized part of the tunnel is loaded. Responses to a stationary unit (1 N) vertical load are calculated for 11 one-third octave band centre frequencies between 20 and 200 Hz.

The exaggerated instantaneous deformation of the cross-section at x=0, containing the load at 200 Hz, is shown in figure 6 for the two designs. The scaling factor for the displacement of the twin bore is five times that for the single bore. At this frequency, strong vibration occurs at the invert in the twin bore tunnel and bending waves propagate either side of the load along the tunnel wall towards the tunnel crown with a high rate of attenuation. On the twin bore design, a harmonic pattern of sixth order in the circumferential direction can be seen. For the single bore, strong bending vibration is observed on the floor, the middle wall and on the tunnel wall below the floor. More vibration is propagated towards the tunnel crown via the middle wall than the tunnel wall itself. This behaviour is also found for lower frequencies.

Figure 6

Mesh and deformation at 200 Hz. (a) For the twin bore design (b) for the single bore design. Dots indicate node points; load shown by triangle symbol.

The vibration strength at a point on the ground surface is quantified by the magnitude of the ‘pseudo-resultant’ displacement of the three components at that point, Embedded Image. Figure 7a shows the displacements on the ground surface immediately above the tunnel crowns for 20, 40, 80 and 160 Hz. At the first three frequencies, the vibration attenuation in the tunnel direction is stronger for the single bore, and therefore, beyond a certain distance the vibration from the single bore is less than that from the twin bore, despite being greater at x=0. This is because the beam-like modes in the tunnels are significant at these frequencies and the single bore tunnel has greater bending stiffness. However, due to the excitation of higher-order bending waves in the tunnel lining at 160 Hz, the vibration attenuations in the two designs are similar and the vibration of the twin bore is less than that of the single bore at all distances.

Figure 7

Pseudo-resultant displacement amplitude (a) on the ground surface immediately above the tunnel crown, (b) along the y-axis. ---, 20 Hz; ---, 40 Hz; -.-, 80 Hz; ⋯, 160 Hz. The heavy lines are for the twin bore and the light lines for the single bore.

Figure 7b presents the displacements along the y-axis on the ground surface for the same four frequencies. The vibration attenuation in the lateral direction is weaker than in the tunnel direction, especially for high frequencies. In general, the maximum response does not occur immediately above the tunnel crown (note: the crown of the twin bore tunnel is at y=10.75 m while that of the single bore is at y=0 m) and the maximum-response position moves towards the tunnel as frequency increases.

In figure 7b, the scattering effect of the companion tunnel in the twin bore design is significant at 20 and 40 Hz, since the response is not symmetrical about the vertical plane containing the axis of the loaded tunnel and the vertical load. At these frequencies the shear wavelengths, 35 and 15 m, are greater than half the distance between the two tunnels. However, at higher frequencies the shear wavelengths are shorter than half the distance between the two tunnels, and the scattering effect of the unloaded tunnel is not seen.

The frequency responses of the ground surface are shown in figure 8. In this figure, displacements at five points on the y-axis are displayed for all 11 frequencies considered. In both cases, broad peaks are found corresponding to resonances in the tunnel structures. However, a much stronger resonance is observed around 40 to 50 Hz for the twin bore when the observer is approximately above the tunnel. This effect may be produced by the wave reflection from the companion tunnel.

Figure 8

Pseudo-resultant displacement amplitude along the y-axis on the ground surface (a) for the twin bore, (b) for the single bore. ---, y=0 m; Embedded Image, y=5 m; ---, y=10 m; ⋯, y=20 m; -.-, y=40 m.

The vibration level is defined as VL=20 log(D/D0) (dB), where D is the pseudo-resultant displacement amplitude and D0 is a reference value. Figure 9a shows the vibration level differences (the level of the twin bore minus that of the single bore) between the two designs for five points immediately above the tunnel crown on the ground surface for the chosen frequencies. Figure 9b shows the vibration level differences for five points along the y-axis on the ground surface. These two figures clearly indicate that, although the loads in these two designs are at the same depth from the ground surface, the twin bore design produces substantially lower levels of vibration than the single bore design. The disadvantage of the single bore in terms of ground vibration may be due to its greater radius, the offset loading (causing the bending of the middle wall which transmits vibration to the tunnel crown) and the shorter distance (3 m) from the tunnel crown to the ground surface.

Figure 9

Vibration level differences (a) on the ground surface immediately above the tunnel crown, (b) along the y-axis on the ground surface. Distances from plan position of the load. ---, 0 m; Embedded Image, 5 m; ---, 10 m; ⋯, 20 m; -.-, 40 m.

6. Moving-load example

The wave speeds in the ground at the depth of the tunnel in the above case are higher than any realistic speed of a railway vehicle. Without recourse to unrealistic speeds, no significant effect can be shown using the method's capability to model the response to a moving harmonic load. A simple example of a load travelling along a track on an embankment on soft alluvial soil is therefore presented. Figure 10 shows the modelled cross-section of the track and ground. The ground is a site at which the properties have been measured (Peplow et al. 1999). It consists of a 2 m layer of weathered alluvial soil with shear wave speed of 72 m s−1, dilatational wave speed of 300 m s−1, a density of 1500 kg m−3 and a damping loss factor of 0.1. A factor such as this, higher than would be expected for soil at depth, has been found to be appropriate for surface vibration propagation (Peplow et al. 1999; Lai et al. 2000). The weathered layer overlies a half-space of stiffer material with shear wave speed of 250 m s−1, dilatational wave speed of 1200 m s−1, a density of 2000 kg m−3 and a damping loss factor of 0.04. The embankment, modelled with FEs, has been attributed the same properties as the substratum. For simplicity the track has been modelled within the profile of the sleepers, using properties that reflect the bending stiffness and mass per unit length of a track panel with UIC rail and typical 300 kg sleepers every 0.6 m.

Figure 10

The modelled half cross-section of the surface track (model is symmetric about the left-hand vertical edge).

Figure 11 shows the amplified response of the ground surface to a 32 Hz load for load speeds of 0 and 70 m s−1 (i.e. close to the shear wave speed in the upper soil layer). The response is plotted at anti-phase with respect to the load which results in the response being shown as upwards at the load position in figure 11a. The load is at the centre position of the track moving towards the right in figure 11b. The surface is modelled to a distance y=30 m and x=±31.25 m. The two responses are scaled by the same factor. Figure 11a shows the waves propagating equally in both directions in the embankment and radially from the embankment into the soil. The response on the soil surface shows an interference pattern because the embankment forms a distributed source. For the moving source, figure 11b shows the response of the embankment distorted asymmetrically and the waves propagating at an angle behind the load. A strict ‘Mach cone’ effect is not seen because of the multiple phase speeds in the ground and the extended nature of the embankment source.

Figure 11

Exaggerated deformation of the ground surface profile at 32 Hz (a) for a stationary load and (b) for a load moving at 70 m s−1.

7. Conclusions

A model using the wavenumber finite and BEMs has been developed for the prediction of responses of structures which are homogeneous in one direction and subject to harmonic loads moving uniformly in that direction. The cross-section of the subdomains may be either open or closed, and may be arbitrarily shaped. Thus, this model is useful for the analysis of many structure–soil interaction problems involving stationary or moving harmonic excitations, such as ground vibration generated by underground trains. Illustrative results from this model for two tunnel designs proposed for an underground railway project show tunnel design significantly affects the ground vibration level. Under the assumptions made, a single bore tunnel large enough to carry two tracks produces substantially higher levels of vibration at the ground surface than two single-track tunnels. The usefulness of the model in showing the effects of a load movement close to the wave speeds in the ground has been illustrated with a simple surface-propagation example.


This work was supported by the Engineering and Physical Science Research Council of the UK under research grant GR/R67309/01.


  • Present address: Holset Engineering Co. Ltd, St Andrews Road, Huddersfield HD1 6RA, UK.

    • Received July 2, 2004.
    • Accepted January 9, 2005.


View Abstract