Modelling of acoustic and elastic wave propagation from underground structures using a 2.5D BEM-FEM approach

Modelling of acoustic and elastic wave propagation from underground structures using a 2.5D BEM-FEM approach A. Romero1, P. Galv´ın1, J. Ant´onio2, J. Dom´ınguez1, A. Tadeu2 1Escuela T´ecnica Superior de Ingenier´ıa, Universidad de Sevilla, Camino de los Descubrimientos, 41092 Sevilla, Spain 2Department of Civil Engineering, University of Coimbra, P´olo II, Rua Lu´ıs Reis Santos, 3030-788 Coimbra, Portugal Abstract This paper presents a numerical method based on a two-and-a-half dimensional (2.5D) boundary elementfinite element (BEM-FEM) coupled formulation to study noise and vibration from underground structures. The proposed model properly represents the soil-structure interaction problem and the radiated noise and vibration. The soil is modelled with the boundary element method, and the Green’s function for a fluid-solid formation is taken as the fundamental solution to represent a solid half-space flattened by a fluid medium, which represents the soil and the air above the ground surface. The finite element method is used to represent structures and enclosed air volumes. The problem representation is limited to a soil-structure interface and the ground surface does not need to be discretised. Radiated noise and vibration are determined after the soil-structure interaction problem has been solved. We verify the proposed method by comparing the solution with an analytical solution for the wave propagation in a fluid-solid medium. Three examples are given to illustrate the noise and vibration radiated by tunnels. The results show that the soil-structure interaction influences the sound pressure field above the ground surface. Keywords: 2.5D BEM-FEM, noise and vibration, tunnels, soil-structure interaction, fluid-structure interaction 1. Introduction In their efforts to achieve sustainable growth many countries have implemented an efficient transportation model to separate economic growth from the use of resources, with lower use of fuel and carbon based systems. Several areas have therefore been modernised, including connected urban and interurban train, metro, and tram systems. However, noise and vibration are among the most common environmental impacts of mass transit solutions, which stem from different sources [1]. Transportation noise can indeed cause annoyance by disrupting sleep, interfering with communication, adversely affecting health and even adversely affecting academic performance. Some studies on general annoyance and activity disturbances show that Email address: [email protected] (A. Romero1) Preprint submitted to Engineering Analysis with Boundary Elements December 17, 2016 train noise is more of a nuisance when there is simultaneous exposure to vibration. This is because it is difficult to distinguish between noise and vibration and this exacerbates the annoyance from noise [2, 3]. Because transport noise pollution affects large numbers of people, railway infrastructure managers require an assessment of the disturbances caused by transport operation and their mitigation. Ground-borne vibration and radiated noise from underground traffic is another major environmental concern in urban areas. The vibration propagates through the underground structure and the surrounding soil and it is perceived directly and sensed indirectly as radiated noise. Technically-supported decisions should be backed up by the accurate numerical modelling of acoustic and elastic waves. The BEM is well suited to represent unbounded problems and the discretization is limited to the boundary. The influence matrices of the BEM are fully-populated and its computation would involve a high computational effort when the problem size is large. Some authors have proposed enhancements based on the fast multipole expansion method (FMM) [4, 5] and on the accelerated BEM by the precorrected fast Fourier transform method (pFFT) [6]. The accelerated BEM based on hierarchical matrices Hand H2 allows the matrix-vector multiplication with almost linear complexity [7, 8]. Brunner et al. [9] compared the BEM based either on the FFM and the hierarchical matrices for the Helmoltz problem, concluding that the computation effort for matrix-vector products was significantly lower when using the hierarchical matrices. Messner and Schanz [10] used H-matrices in the development of an accelerated time-domain BEM for elasticity. These authors separated the far-field from the near-field by means of the hierarchical matrices. They improved the compression of the BEM matrices by separating into a sum of Hand H2matrices. The accelerated methodologies reduce the computational effort in order to represent three-dimensional problems but the 2.5D formulation is preferred to study quite long invariant geometries. For this, several numerical models based on two-and-a-half dimensional (2.5D) formulations have been proposed to account for the longitudinally invariant geometry of structures such as tunnels. These formulations compute the three-dimensional (3D) wave field from two-dimensional (2D) problems with different wavenumbers [11] to avoid the computationally expensive disadvantage of 3D formulations. Forrest and Hunt [12] are among the authors who have proposed a semi-analytical solution for the timeharmonic displacements of a tunnel modelled as an infinitely long, thin cylindrical shell, where the soil was represented as a homogeneous full-space. The coupled problem was solved in the frequency domain by Fourier decomposition into ring modes circumferentially and a Fourier transform into the wavenumber domain longitudinally. Fran¸cois et al. [13] presented a 2.5D coupled BEM-FEM methodology to compute the dynamic interaction between a layered soil and structures with longitudinally invariant geometry. This formulation uses a 2.5D Green’s function for a layered half-space, thus there is no need to discretise the free surface and the layer interfaces. A regularised boundary integral equation is derived to avoid the evaluation of singular traction integrals. Later, Galv´ın et al. [14] used that formulation to predict railway induced vibrations in a tunnel embedded in a layered half-space. Moreover, Gupta et al. [15] proposed a coupled 2 periodic BEM-FEM approach, where a boundary element method is used for the soil and a finite element method for the tunnel. The tunnel periodicity is handled by the Floquet transformation to formulate the tunnel-soil interaction problem in the frequency-wavenumber domain, and to compute the wave propagation field in the soil. Furthermore, a windowing technique was proposed by Coulier et al. [16] to represent finite domain dimensions in a 2.5D formulation to study SSI problems. The proposed method allowed the study of finite insulation trenches within an invariant soil geometry. Many works suggest radiated noise should be computed assuming that the sound receiver has no effect on the vibration generation mechanism. Therefore, this simplification allows the decoupling of the incident wave field and the radiated noise. Following this procedure, Nagy et al. [17] proposed a Rayleigh integral-based method combined with the FEM to predict radiated noise in buildings from underground railway traffic. This model was experimentally validated in a building close to a railway line in Paris and maximum sound pressure level was found to be 75 dB. Fiala et al. [18] also suggested a numerical model to study vibrations and radiated noise caused by underground railways. This model solves the soil-structure interaction problem with a BEM-FEM coupled formulation using the Green’s function for the layered half-space. The acoustic radiation problem is then computed assuming weak coupling between structural and acoustic waves. An acoustic spectral finite element method is used to predict the radiated noise. Tadeu et al. [19] subsequently modelled the acoustic attenuation provided by a barrier in an underground train station using a comprehensive coupled formulation based on the BEM and the method of fundamental solutions (MFS). In this case, the proposed method can compute the noise without the former simplifications. Recent publications have produced computationally feasible numerical models which can represent the ground-borne vibration at specific sites. Romero et al. [20] presented a 2.5D BEM-FEM model to study noise and vibration within a tunnel embedded in an unbounded solid due to a moving load. They used the BEM for the full-space fundamental solution in elastodynamic and the FEM in fluid-acoustics and elastodynamics. The results showed that the tunnel displacement and the air pressure inside the tunnel increase with the load speed according to the regime defined by the wave propagation velocity in each medium. The same behaviour was observed in the unbounded solid. However, to the best of our knowledge the radiated noise and vibration at the soil surface from underground structures has not been studied using a coupled formulation. This work describes a novel numerical model that can handle the above mentioned problem. The model takes a domain decomposition approach to study fluid acoustics and solid scattering waves in a half-space medium with an innovative BEM technique. The BEM formulation considers, as fundamental solution, a problem where a formation composed of an elastic solid medium bounded by an acoustic fluid medium that can be the soil as well the air above the ground surface. Therefore, the discretisation of the fluid-solid interface is not needed and only soil-structure interfaces are modelled. The proposed methodology is used to study the acoustic and elastic scattered wave field inside underground structures, within the soil, and above the ground surface. The paper is organised as follows. First, the numerical model is presented, the 2.5D 3 BEM-FEM fluid-acoustics formulation is described, and the Green’s function for a solid half-space bounded by a fluid formation is discussed. The proposed method is then verified with a benchmark problem that studies the wave propagation in fluid-solid formations. This section also analyses the convergence of the longitudinally spatial Fourier transform. The proposed method is then used to study noise and vibration in three examples related to the wave field radiated by a single tunnel, a twin tunnel, and an underground station with a sound insulation screen. Finally, the conclusion section summarises the main contributions of this work. 2. Numerical model The numerical model is based on a coupled BEM-FEM formulation (Figure 1). The underground structure, composed by solid (Ωs) and fluid (Ωf) enclosures, is represented by the FEM. The interface between both enclosures is denoted as Γq. The soil domain (Ωs∞) and the air above the ground surface (Ωf∞) are described by the BEM. The Green’s function for a fluid-solid formation presented by Tadeu and Ant´onio [21] is used as fundamental solution in the BEM. Therefore, the boundary element discretisation is limited to the limiting interface between the soil and the structure (Γo). The coupled BEM-FEM formulation is addressed by imposing proper conditions at both interfaces, the fluid-structure (Γq) and the soil-structure (Γo) interfaces. The radiated wave field is computed once the displacement and traction solutions at Γoare known. Ωf∞ Ωs∞ Γo Γq Ωf Ωs y x z Figure 1: Problem definition (figure not scaled). Solid (Ωs) and fluid (Ωf) enclosures, soil domain (Ωs∞), air above the ground surface (Ωf∞) and the fluid-structure (Γq) and the soil-structure (Γo) interfaces. The 2.5D formulation computes the 3D solution in the frequency domain, assuming that the problem is invariant in the longitudinal direction z, as the superposition of 2D problems with a different longitudinal wavenumber (kz) in the zdirection. An inverse Fourier transform is used to compute the 3D solution at a 4 point x(x, y, z): a(x, ω) = Z+∞ −∞ ba(b x, kz, ω)e−ikzzdkz(1) where a(x, ω) is the unknown variable in the frequency domain (e.g. displacement, traction or pressure), ba(b x, kz, ω) is its representation in the wavenumber-frequency domain at b x=x(x, y, 0), ωis the corresponding angular frequency, and i=√−1. The Equation (1) converges and can be approximated by a discrete Fourier transform over kz. The wavenumber sampling ∆kzdefines the maximum distance zwhere the solution is accurately computed as max(z) = 2π/∆kz/2. In this work, the wavenumber sampling used in this work is uniform spaced. Moreover, the problem solution is obtained at equally spaced points with ∆z= 2π/ max(kz). The computational effort of the spatial Fourier transform depends on these parameters. Next, the numerical methodology is discussed. First, the BEM formulation and its fundamental solution are presented. Then, the FEM formulation for elastodynamic and fluid-acoustics are briefly described. Finally, the coupled BEM-FEM model is obtained. 2.1. Boundary element formulation The boundary element formulation we present considers an arbitrary boundary inside a solid half-space. This solid half-space is bounded by a fluid medium. The solid and the fluid half-spaces are used to represent the soil domain (Ωs∞) and the air above the ground surface (Ωf∞), respectively. The arbitrary enclosure accounts for the soil-structure interface (Γo). The integral representation of displacement in the frequency domain for a point xilocated at the arbitrary boundary Γocan be written as [22]: ci lkui k(xi, ω) = ZΓo u∗ lk(x, ω;xi)ti k(x, ω)dΓ−ZΓo t∗ lk(x, ω;xi)ui k(x, ω)dΓ(2) where ui kand ti kare the component kof displacements and tractions at the interface Γo, respectively. u∗ lk and t∗ lk are the fluid-solid domain fundamental solution for displacements and tractions, respectively, at point x due to a point load at xiacting in the ldirection at the solid [21]. The integral-free term ci lk depends only on the boundary geometry at point xi[22]. Assuming that the boundary Γois invariant in the longitudinal direction z, Equation (2) is expressed in terms of integrals over the cross section of the boundary Γo, denoted as Σo: ci lkui k(xi, ω) = ZΣoZ+∞ −∞ u∗ lk(x, ω;xi)ti k(x, ω)dΣdz −ZΣoZ+∞ −∞ t∗ lk(x, ω;xi)ui k(x, ω)dΣdz (3) Equation (3) is then transformed to the wavenumber domain as (Equation 1): 5 ci lk bui k(b xi, kz, ω) = ZΣobu∗ lk(b x, kz, ω;b xi)bti k(b x, kz, ω)dΣ−ZΣobt∗ lk(b x, kz, ω;b xi)bui k(b x, kz, ω)dΣ(4) The boundary cross section Σois discretised into elements and nodes, leading to a boundary approximation of the displacement and traction using the interpolation shape function N. Then, Equation (4) is written as: ci lk bui k= Q X j=1 "(ZΣjbu∗ lkNjdΣ)btj k−(ZΣjbt∗ lkNjdΣ)buj k#= Q X j=1 hb Gij lkbtj k−b Hij lk buj ki(5) where Qis the number of nodes at Σo, and Σjstands for the elements which contain the node j. The system of equations for all the boundary nodes in the wavenumber-frequency domain becomes: b H(kz, ω)b uo(b x, kz, ω) = b G(kz, ω)bto(b x, kz, ω) (6) where b Gand b Hare the fully nonsymmetrical boundary element system matrices, and b uoand btoare displacements and tractions at the interface Γo, repectively. In this work, constant boundary elements are used. The boundary integrals are computed using an eight point Gauss-Legendre quadrature whenever the collocation point does not belong to the integration element. Closed forms are used to solve singular integrals when the collocation point b xibelongs to the integration element. The exact expressions for the singular integrals given in References [23, 24] are used. The analytical solution for the steady state response of a formation composed by a homogeneous solid half-space bounded by a fluid medium subjected to a spatially harmonic load can be found in Reference [21]. This solution relates displacement at point x(x, y, z) due to a dynamic source placed at a different point xi(xi, yi, zi) in the fluid-solid formation (Figure 2). The boundary element collocation method described by Equation (5) requires evaluation of the fundamental solution at each Gaussian integration point for sources placed at every nodal collocation point. This procedure could be highly expensive, depending on the size of the soil-structure interface discretisation. In this section, the Green’s function presented in Reference [21] is modified in a feasible form to carry out the BEM collocation methodology. The Green’s function (G) given by a spatially sinusoidal harmonic load in the zdirection in the frequency domain is expressed as the superposition of 2D problems with the longitudinal wavenumber kzvarying in the zdirection: G(x, ω;xi) = Z+∞ −∞ b G(b x, kz, ω;b xi)e−ikz(z−zi)dkz(7) where b Grepresents the Green’s function in the wavenumber-frequency domain. The Green’s function for the fluid-solid formation is derived by imposing compatibility and equilibrium of normal displacements and stresses, and assuming null shear stresses, at the interface between the solid 6 y x z Ωf∞ Ωs∞ φi f φfφi s, ψi x, ψi y, ψi z φs, ψx, ψy, ψz xix xix Figure 2: Solid formation bounded by a flat fluid medium (figure not scaled). Soil domain (Ωs∞), air above the ground surface (Ωf∞), source point (xi(xi, yi, zi)), observation point (x(x, y, z)), solid compressional (φs,φi s) and rotational (ψx,ψy,ψz,ψi x, ψi y,ψi z) potentials and fluid dilatational potential (φf,φi f). and the fluid medium. The function can be expressed as the superposition of: (i) source terms equal to those in the full-space ( b Gfull), and (ii) surface terms ( b Gsurf ) needed to represent the indirect wave field scattered by the fluid-solid interface [21]: b G(b x, kz, ω;b xi) = b Gfull(b x, kz, ω;b xi) + b Gsurf (b x, kz, ω;b xi) (8) The expressions for the Green’s function for the full-space ( b Gfull) are defined explicitly in Reference [25]. The source terms are obtained from compressional (φs) and rotational (ψx,ψy,ψz) potentials for a solid full-space. The fluid dilatational potential (φf) for a harmonic pressure load acting on an unbounded fluid medium defines the Green’s function for the direct pressure field. The meaning of these potentials is well established in Reference [21]. The surface terms can be expressed in the same way as those of the source terms using the potentials φi s,ψi x,ψi y,ψi z, and φi f. Each potential amplitude can be computed imposing the four stated boundary conditions at the fluid-solid interface (y= 0):                    σs yx =σf yx = 0 σs yz =σf yz = 0 σs yy =σf yy us y=uf y (9) where superscripts sand frelate to variables belonging to the solid and fluid subdomains, respectively. 7 The surface terms are given by a spatial Fourier transform in the xdirection for the wavenumber kx, which can be approximated by a finite sum of terms: b Gsurf (b x, kz, ω;b xi) = Z+∞ −∞ e Gsurf (kx, y, kz, ω;yi)e−ikx(x−xi)dkx(10) Following the approach given in Equation (10), the expressions for solid displacements due to a load acting in the xdirection are defined as [21]:              e Gsurf xx =EaAx−ik2 x νs x Eb+−iγxCx−ikz γx BxEc e Gsurf yx =Ea(−ikxAxEb+ikxCxEc) e Gsurf zx =Ea−ikzkx νs x AxEb+ikzkx γx BxEc(11) where Ea= 1/2ρsω2,Eb=e−iνs xy,Ec=e−iγxy,νs x=qk2 p−k2 z−k2 xwith (Im(νs x)≤0), and γx= pk2 s−k2 z−k2 xwith (Im(γx)≤0). The wavenumbers related to P and S waves are kp=ω/cpand ks=ω/cs, respectively, ρsis the solid density, and cpand csare the propagation velocity of P and S waves in the solid, respectively. Ax,Bx, and Cxare unknown coefficients to be determined from the conditions defined by Equation (9). The imposition of the four boundary conditions leads to the following two systems of four equations [21]:           −2k2 x−k2 zk2 x−γx20 −2 1 1 0 −k2 s−2ν2 zx νs x 0 2γx−i2ρsω2 νf xµs −i0i−2ρsω2 k2 fλf                   Ax Bx Cx Dx         =          −2k2 x−k2 s+ 2k2 x −2 2 k2 s+ 2ν2 zx νs x−2γx −i i          (12) where µsis the solid second Lam´e constant, λfis the fluid Lam´e constant, νf x=qk2 f−k2 z−k2 xwith (Im(νf x)≤0), νzx =p−k2 z−k2 x, and the wavenumber related to the fluid pressure wave is kf=ω/cf, with cfbeing the fluid compressional wave velocity. After the solution vectors Ax,Bx,Cxand Dxare computed, the amplitudes of the potentials are given by the following expressions:            Ax=Ei bAx(1) + Ei cAx(2) Bx=Ei bBx(1) + Ei cBx(2) Cx=Ei bCx(1) + Ei cCx(2) (13) with Ei b=e−iνs xyiand Ei c=e−iγxyi. The fluid pressure is computed from Dx[21]. 8 In the same way, the fundamental solution for a load acting in the ydirection is expressed through the following equations:              e Gsurf xy =Ea[−iAykxEb+iBykxEc] e Gsurf yy =Ea−iνs xAyEb+−ik2 x γx By+−ik2 z γx CyEc e Gsurf zy =Ea[−iAykzEb+iCykzEc] (14) where the amplitudes of potentials are defined by:             −2νs x−k2 x γx +γx−k2 z γx 0 −2νs x−k2 x γx −k2 z γx +γx0 −k2 s−2ν2 zx −2k2 x−2k2 z−i2ρsω2 νf xµs −iνs x−ik2 x γx −ik2 z γx −i2ρsω2 k2 fλf                     Ay By Cy Dy         =            2νs x−ν2 zx γx +γx 2νs x−ν2 zx γx +γx −k2 s−2ν2 zx2ν2 zx iνs xik2 x+ik2 z γx            (15) and:            Ay=Ei bAy(1) + Ei cAy(2) By=Ei bBy(1) + Ei cBy(2) Cy=Ei bCy(1) + Ei cCy(2) (16) Furthermore, solid displacements for a load acting in the zdirection are obtained as follows:              e Gsurf xz =Ea−ikzkx νs x AzEb+ikzkx γx CzEc e Gsurf yz =Ea[−ikzAzEb+iBzkzEc] e Gsurf zz =Ea−ik2 z νs x AzEb+−ik2 x γx Cz−iγxBzEc(17) The amplitudes of the potentials are provided by the following system of equations:           −2 1 1 0 −2k2 zk2 z−γx2−k2 x0 −k2 s+ 2ν2 zx νs x 2γx0−i2ρsω2 νf xkzµs −ikzikz0−2ρsω2 k2 fλf                   Az Bz Cz Dz         =          −2 2 −2k2 zk2 z−γx−k2 x k2 s+ 2ν2 zx νs x−2γx −ikzikz          (18) thus,            Az=Ei bAz(1) + Ei cAz(2) Bz=Ei bBz(1) + Ei cBz(2) Cz=Ei bCz(1) + Ei cCz(2) (19) 9 kzcs/ω. Superimposed on this figure we can see the dispersion curve kR=cs/cRof the Rayleigh waves with propagation velocity cR= 2442 m/s, and the dispersion curve kF=cs/cfof the dilatational pressure waves in the fluid. The dispersion curve kC=cs/C of waves with phase velocity Cis also represented. The phase wave velocity was estimated as: C(ω) = ω pk∗2 x(ω) + k∗2 z(ω)(38) where k∗ xand k∗ zdenotes the wavenumbers where e Gsurf yy (kx, y, kz, ω;b x0) takes its maximum value. (a) (b) (c) (d) 0 10.50.25 0.75 Figure 5: Real part of the Green’s function for displacements (a) e Gsurf xy (×1014 m), (b) e Gsurf yy (×1014 m), and (c) e Gsurf zy (×1014 m) at receiver R1, and pressure (d) eσfy (×105Pa) at receiver R2. Dimensionless wavenumber kR(solid line), kC (dashed line), and kF(dotted line) are also represented. The aforementioned dispersion curves define circumferential geometric regions where surface terms reach local maxima. Afterwards the Green’s function varies smoothly and decays to zero value. In this work, the surface terms e Gsurf were computed using 500 logarithmic spaced points for kx, from 10−6to 106, to ensure the convergence of Equation (10). Figure 6 compares the analytical solution and the numerical results for solid displacement at receiver R1, and the pressure field at receiver R2. A constant wavenumber kz= 0.4 rad/m was assumed for the verification 16 of the proposed formulation. Numerical results are in good agreement with the reference solution [21]. 0 40 80 120 160 200 240 280 320 −3 −1.5 0 1.5 3x 10−12 Frequency [Hz] Displacement [m] (a) 0 40 80 120 160 200 240 280 320 −2 −1 0 1 2x 10−11 Frequency [Hz] Displacement [m] (b) 0 40 80 120 160 200 240 280 320 −4 −2 0 2 4x 10−12 Frequency [Hz] Displacement [m] (c) 0 40 80 120 160 200 240 280 320 −0.04 −0.02 0 0.02 0.04 Frequency [Hz] Pressure [Pa] (d) Figure 6: Real (black) and imaginary (grey) parts of the analytical solution [21] (solid lines) and numerical results (circles) for (a) horizontal, (b) vertical, and (c) longitudinal displacements at receiver R1, and (d) pressure at receiver R2, due to a vertical load acting at So. The 3D problem solution is given by Equation (1) as the superposition of 2D problems with different wavenumbers kz. This equation converges and can be solved by a finite sum of terms. In order to study the convergence, Figure 7 shows the problem solution at R1and R2for a range of frequencies and wavenumbers kz, and the previous dispersion curves were also represented. The phase wave kCand the Rayleigh wave kRdispersion curves intersect each other at the frequency fco = 118 Hz. Local maxima were found around kCat lower frequencies, after which the solution decays for wavenumbers higher than kR. Additionally, the fluid pressure is mainly found around the region defined by kRand kF, and becomes significant from the frequency. Thus, the convergence of Equation (1) would be ensured for wavenumbers kzhigher than kFfrom the frequency fco. However, the computation for lower frequencies should be extended for a few wavenumbers until the solution decays to zero. 17 (a) (b) (c) (d) 0 10.50.25 0.75 Figure 7: Amplitude of (a) horizontal (×1012 m), (b) vertical (×1011 m), and (c) longitudinal (×1012 m) displacements at receiver R1, and (d) pressure field (×102Pa) at receiver R2, due to a vertical load acting at So. Dimensionless wavenumber kR(solid line), kC(dashed line), and kF(dotted line) are also represented. 4. Applications This section explores ground vibrations and radiated noise in three applications. First, the wave field radiated by a single tunnel due to harmonic loads is evaluated. Next, the scattered wave field caused by the insertion of a second tunnel is studied in a twin tunnel system. Finally, changes in sound pressure at the free-field produced by the insertion of mitigation systems in the underground structure are assessed. In this last example, an underground railway station with an acoustic barrier is considered. 4.1. Radiated noise and vibration generated by a tunnel This example analyses the noise and vibration field radiated from a tunnel embedded in a homogeneous half-space (Figure 8). Soil displacements and the sound pressure level inside the tunnel and above the ground surface due to a harmonic load P= 2πN are studied for a maximum frequency range of fmax = 250 Hz. The soil response of this problem has previously been studied by Gupta et al. [31]. This reference has been given to verify the conclusions obtained with the proposed methodology. 18 y x z Ωf∞ Ωs∞ P Γo Γq Ωs Ωf + d Figure 8: Definition of tunnel and subdomains (figure not scaled). Soil domain (Ωs∞), air above the ground surface (Ωf∞), solid (Ωs) and fluid (Ωf) enclosures, the fluid-structure (Γq) and the soil-structure (Γo) interfaces, tunnel depth (d) and acting force (P). The tunnel is at depth d= 20 m, has radius r= 3.0 m and wall thickness t= 0.3 m. The lining has concrete properties with Young’s modulus E= 35 ×109N/m2, Poisson’s ratio ν= 0.25, and density ρ= 2500 kg/m3. The soil has a P-wave propagation velocity cp= 500 m/s, an S-wave propagation velocity cs= 250 m/s, and a density ρs= 1750 kg/m3. The sound propagation velocity is cf= 340 m/s, and the air density takes a value of ρf= 1.22 kg/m3. The analysis was carried out in the frequency-wavenumber domain for frequencies ranging from 2 Hz to 250 Hz, with a frequency step of 2 Hz. The longitudinal wavenumber varies from kz= 0 to kz= 3×kF= 4.08, which ensure the radiated field is computed accurately at receivers far from the application load. The wavenumber step was set so as to adequately represent n= 100 acoustic wavelengths, ∆kz= 1/n. The tunnel lining is represented with 216 shell finite elements (Ωs), and the discretisation matches with the boundary element mesh used for the soil-structure interface (Γo), as well as the fluid-structure interface (Γq). The 2.5D shell element is a two-nodes line element with four degree of freedom at each node (three displacements and one rotation around the zaxis). The air volume inside the tunnel (Ωf) is represented with 9562 fluid finite elements [20], allowing the representation of wavelengths with at least 6 elements. Figure 9 shows vertical displacements and sound pressure levels at three observation points on the ground surface (y= 0 m), at different distances from the tunnel axis: x= 0 m, x= 4 m, and x= 16 m. The frequency content has an undulating behaviour because different waves are propagating through the fluid-solid formation, and the frequency step between subsequent undulations becomes lower as the distance of the observation point increases, but remains uniform in the frequency range. The Green’s function for the fluid-solid formation is also represented in this figure to show that differences between the field radiated by the tunnel and the fluid-solid formation behaviour are due to the waves scattered by the structure. The 19 mismatch becomes more significant as fluid-solid wavelengths are smaller than the tunnel size. These results are consistent with those reported in Reference [31] for the soil response. Frequency [Hz] 0 50 100 150 200 250 Displacement [dB ref 1m/N] -280 -260 -240 -220 -200 -180 (a) Frequency [Hz] 0 50 100 150 200 250 Displacement [dB ref 1m/N] -280 -260 -240 -220 -200 -180 (b) Frequency [Hz] 0 50 100 150 200 250 Displacement [dB ref 1m/N] -280 -260 -240 -220 -200 -180 (c) Frequency [Hz] 0 50 100 150 200 250 Sound pressure level [dB ref 20 µPa/N] -80 -60 -40 -20 0 (d) Frequency [Hz] 0 50 100 150 200 250 Sound pressure level [dB ref 20 µPa/N] -80 -60 -40 -20 0 (e) Frequency [Hz] 0 50 100 150 200 250 Sound pressure level [dB ref 20 µPa/N] -80 -60 -40 -20 0 (f) Figure 9: (a-c) Vertical displacement (black line), and (d-f) sound pressure level (black line), at points on the ground surface with (a,d) x= 0 m, (b,e) x= 4 m, and (c,f) x= 16 m due to a load acting at the tunnel invert. The analytical solution for a homogeneous fluid-solid formation [21] is also represented (grey line). Vertical soil displacement and sound pressure level distribution, both inside the tunnel and above the ground surface, are represented in Figure 10 for frequencies of 10 Hz and 80 Hz. Results were computed at a grid of 7141 receivers equally spaced, and the 3D solution was computed afterwards. Maximum displacements were found around the tunnel invert according to the soil wavelength. A shadowed part is observed above the tunnel where the response is lower. Moreover, the sound pressure distribution is defined by the fluidmomentum relation, which links fluid pressure with soil displacement. Therefore, radiated waves at the free field show a maximum above the tunnel, and the maximum inside the tunnel is found around the lining. Fluid wavelengths are longer than those in the solid, according to the wave propagation velocity in each medium. 4.2. Scattered wave field by a twin tunnel system This example examines the wave field scattered by two tunnels embedded in a homogeneous half-space (Figure 11). One of these tunnels is subject to a dynamic force P= 2πN applied to the tunnel invert. The importance of the tunnel’s interactions is analysed through the insertion gain that results from the existence 20 (a) (b) −110 −0.50.5 Figure 10: Real part of tunnel and soil vertical displacements (×107m), and sound pressure inside the tunnel (×102Pa at 10 Hz and ×101Pa at 80 Hz) and above ground surface (×103Pa), produced by a point load of (a) 10 Hz and (b) 80 Hz acting on the tunnel invert. of a second tunnel. Kuo et al. [32] studied the effect of a twin tunnel system on the wave propagation in an unbounded soil domain. The same system is analysed in this example. The radiated noise is also studied. y x z Ωf∞ Ωs∞ P c ΩfΩf Γo Γq Ωs d + + Figure 11: Diagram of twin tunnel (figure not scaled). Soil domain (Ωs∞), air above the ground surface (Ωf∞), solid (Ωs) and fluid (Ωf) enclosures, the fluid-structure (Γq) and the soil-structure (Γo) interfaces, tunnel depth (d) and acting force (P). Both tunnels are located at depth d= 20 m and they are separated by a distance of c= 10 m, they have radius r= 3.0 m and the wall thickness is t= 0.3 m. The lining has the properties of concrete, with 21 Young’s modulus E= 35 ×109N/m2, Poisson’s ratio ν= 0.25, and density ρ= 2500 kg/m3. The soil has a P-wave propagation velocity cp= 500 m/s, an S-wave propagation velocity cs= 250 m/s, and a density ρs= 1750 kg/m3. The sound propagation velocity is cf= 340 m/s, and the air density takes a value of ρf= 1.22 kg/m3. The response of the twin-tunnel system is computed using the proposed methodology. Each tunnel is modelled with 56 shell elements and the volume of air within the tunnels is discretised with 588 fluid elements. The boundary element mesh matches the finite element discretisation. The solution is computed for a vertical load applied on the left-hand tunnel invert (Figure 11) acting at 60 Hz. The wavenumber kz range was set as in the previous example. Figure 12.(a) shows the vertical soil displacements and the sound pressure inside tunnels. The soil response exhibits a non-symmetric distribution around the vertical axis of the loaded tunnel due to the wave field scattered by the right-hand tunnel. Maximum displacements occurred around the source point. The sound pressure pattern of both tunnels is defined by the soil displacement around them. The sound pressure reaches its highest value in the loaded tunnel. (a) (b) Figure 12: (a) Real part of vertical soil displacement (×107m) and radiated noise inside the tunnels (×101Pa), due to a harmonic load acting on the tunnel invert at f= 60 Hz, and (b) the insertion gain from the response of a twin-tunnel system over that of a single-tunnel. The effect of waves scattered by the right-hand tunnel is investigated in terms of the insertion loss in the single-tunnel and twin-tunnel systems. The second tunnel changes both the soil displacement distribution and the sound pressure inside the loaded tunnel, and this becomes more important at the right-hand tunnel, as can be seen in Figure 12.(b). These conclusions agree with those presented by Kuo et al. [32] for the ground response. In fact, the second tunnel modifies the sound pressure level distribution at the ground 22 surface (Figure 13). The sound pressure field is distributed around the source point and comprises a circular region about 20 m. The insertion loss also shows major differences between the tunnel systems. (a) (b) Figure 13: (a) Real part of sound pressure level [dB ref 20µPa] at the ground surface due to a harmonic load acting on the tunnel invert at f= 60 Hz, and (b) the insertion gain [dB] from the response of a twin-tunnel system over that of a single tunnel. Tunnel locations are marked by a dashed line. 4.3. Noise radiated from an underground train station with acoustic insulation The last application examines the noise radiated by an underground railway station (Figure 14) that has an acoustic barrier that separates the two railway tracks. An acoustic point pressure load, located at S0(−2,14) simulates the noise of a passing train. This example was proposed by Tadeu et al. [19] for an unbounded domain to study the sound wave propagation inside the station in 2D cases. We have extended the problem to evaluate the influence of an acoustic barrier on the free-field radiated noise. The station is at a depth d= 12 m and the geometry can be found in Figure 14. The structure has a concrete wall of thickness t= 0.3 m, Young’s modulus E= 50×109N/m2, Poisson’s ratio ν= 0.3, and density ρ= 2500 kg/m3. The soil has a P-wave propagation velocity cp= 500 m/s, an S-wave propagation velocity cs= 250 m/s, and a density ρs= 1750 kg/m3. The air takes sound propagation velocity cf= 340 m/s, and density ρf= 1.22 kg/m3. The acoustic barrier was assumed to be rigid, and 3.5 m long, and it was modelled as a discontinuity for the air volume Ωf. There is a gap of 1 m exists between the station invert and the barrier. The vibration and radiated noise produced by a moving pressure load Pacting at point Soare analysed next. This load reproduces a moving source travelling in the zdirection at speed cf=ω/kz. This speed is not a realistic traveling train speed but its interest is related with the wave propagation phenomenum in mediums with different wave propagation velocities. The source is assumed to be represented as a Ricker wavelet with a characteristic frequency of 145 Hz and amplitude 2πPa [19]. The solution was computed for a frequency range from 2 Hz to 512 Hz, with a frequency step of 2 Hz. Complex frequencies with an imaginary part of the form Im(ω) = −0.7∆ωwere used to reduce the contribution of the virtual source to the response and to prevent the occurrence of aliasing phenomena. Each frequency step was solved for a wavenumber 23 y x z Ωf∞ Ωs∞ Ωf So 4 m r1= 8 m r2= 4 m 2 m d + Γo Γq Ωs Figure 14: Diagram of underground station and definition of subdomains (figure not scaled). Soil domain (Ωs∞), air above the ground surface (Ωf∞), solid (Ωs) and fluid (Ωf) enclosures, the fluid-structure (Γq) and the soil-structure (Γo) interfaces, depth (d) and source point (So). defined by the load speed. The characteristic element size enabled the minimum wavelength to be suitably represented with at least six elements. Figure 15 compares the sound pressure distribution inside the station with and without the insulation barrier at the cross section defined by the coordinate plane y= 12 m. The pressure distribution was represented for the time when the load passes at z= 0 m. Time solutions were obtained by applying an inverse Fourier transform to the frequency response. The sound pressure fields exhibit a typical Mach cone and many differences between the two cases were found due to the presence of the barrier. The acoustic screen induced several wave reflections at the station wall and at the barrier. The free-field radiated noise is due to ground-borne vibration defined by the stated boundary conditions at soil surface (Equation (9)). The sound pressure changes after the insertion of the acoustic barrier, as can be seen in Figure 16, and so does the soil displacement (Figure 17). The barrier redistributes displacements and tractions at the soil-structure interface and this modifies the wave field at the ground surface. The arrival time of S waves gives a spatial delay of 18 m in the displacements and sound pressure field. 5. Conclusions This work has proposed a BEM-FEM formulation to study fluid and solid wave scattering in half-space formations. The method was formulated in 2.5D but is suitable for 3D problems whose material and 24 (a) (b) −550 −2.52.5 Figure 15: Radiated noise inside the station at a plane defined by the coordinate y= 12 m, due to a moving load travelling at v=cf, (a) without acoustic barrier, and (b) with acoustic barrier (black line). (a) (b) −0.01 0.010 −0.005 0.005 Figure 16: Radiated noise at the soil surface due to a moving load travelling at v=cf(a) without acoustic barrier, and (b) with acoustic barrier. geometric properties are homogeneous in one direction. The proposed model has been developed to analyse the acoustic and elastic wave propagation from underground structures. The FEM modelling structures and enclosed fluid subdomains, while the BEM represents the soil-structure interface. The model was verified through a benchmark problem with a known analytical solution, and numerical results were in good agreement with the reference solution. The main practical application of the method is the analysis of the noise and vibration radiated from tunnels, with a comprehensive formulation that predicts both sound pressure and soil displacement. Three examples were given that show the free-field radiated noise dependence on the scattered waves within the 25