Acoustic waves scattered by elastic waveguides using a spectral approach with a 2.5D coupled boundary-finite element method

Acoustic waves scattered by elastic waveguides using a spectral approach with a 2.5D coupled boundary-finite element method F.J. Cruz-Mu˜noza, A. Romeroa,∗, P. Galv´ına, A. Tadeub,c aEscuela T´ecnica Superior de Ingenier´ıa, Universidad de Sevilla, Camino de los Descubrimientos, 41092 Sevilla, Spain bITeCons - Institute for Research and Technological Development in Construction, Energy, Environment and Sustainability, Rua Pedro Hispano s/ n., 3030-289 Coimbra, Portugal cADAI - LAETA, Department of Civil Engineering, University of Coimbra, P´olo II, Rua Lu´ıs Reis Santos, 3030-788 Coimbra, Portugal Abstract This work presents a two-and-a-half dimensional (2.5D) spectral formulation based on the finite element method (FEM) and the boundary element method (BEM) to study wave propagation in acoustic and elastic waveguides. The analysis involved superposing two dimensional (2D) problems with different longitudinal wavenumbers. A spectral finite element (SFEM) is proposed to represent waveguides in solids with arbitrary cross-section. Moreover, the BEM is extended to its spectral formulation (SBEM) to study unbounded fluid media and acoustic enclosures. Both approaches use Lagrange polynomials as element shape functions at the Legendre-Gauss-Lobatto (LGL) points. The fluid and solid subdomains are coupled by applying the appropriate boundary conditions at the limiting interface. The proposed method is verified by means of two benchmark problems: wave propagation in an unbounded acoustic medium and the scattering of waves by an elastic inclusion. The convergence and the computational effort are evaluated for different h−pstrategies. Numerical results show good agreement with the reference solution. Finally, the proposed method is used to study the pressure field generated by an array of elastic fluid-filled scatterers immersed in an acoustic medium. Keywords: SBEM, SFEM, fluid-solid interaction, waveguide, scattered waves, two-and-a-half dimension 1. Introduction Many engineering fields involve time harmonic wave propagation, such as fluid acoustics and solid scattering. Hybrid methods based on the boundary element method (BEM) and the finite element method (FEM) are suitable for studying solid and fluid interaction in unbounded media. Both methods have been used in various works to predict the response of coupled fluid-structure problems. For the low frequency5 range, the standard formulation using linear elements accurately represents the fluid and solid scattering ∗Corresponding author Email address: [email protected] (A. Romero) Preprint submitted to Engineering Analysis with Boundary Elements 29th April 2019 waves. However, at high frequencies, the solution deteriorates due to the so-called pollution effects [1, 2]. Higher-order approximations are needed to obtain accurate results. Different approaches have been proposed in the past to overcome the pollution effect in two and three dimensional problems [3]. The h-refinement method has proved to be useless at high wavenumbers because10 of the necessary computational effort to represent small wavelengths [4]. Many works suggest using highorder elements to improve the convergence rate of standard formulations. The boundary and finite element methods have thus been extended to their spectral formulations [5, 6] approximating the field variables by high-order interpolation shape functions. Such functions are polynomials over non-uniformly spaced nodes to avoid the Runge phenomenon.15 A number of authors have successfully applied spectral finite element methods in different fields. Degrande and De Roeck [7] developed a high-order formulation to study the dynamic response of a porous saturated medium using spectral elements. The authors used the analytical solution of Biot equations in two dimensions as shape functions. Bar-Yoseph et al. [8] proposed a spectral element method to study the transient response of a simply supported Euler-Bernoulli beam subjected to external forced lateral vibrations.20 They examined two Galerkin formulations using Hermitian and Lagrangian polynomials as interpolants for spatial and time discretisations. The space-time SFEM enabled the use of larger time steps while maintaining the accuracy of the solution. Kudela et al. [13] applied a SFEM to study the wave propagation in an isotropic rod and in a Timoshenko beam. The authors proposed this methodology to detect small damage in structures. The computed results were compared with those obtained from a standard FEM and from25 experimental measurements. This comparison highlighted the efficiency of the spectral elements. Zhu et al. [14] proposed the use of Chebyshev spectral elements within an implicit Newmark method to simulate acoustic wave propagation. They found accurate results with uniform nodal element distribution and small time steps. Recently, Romero et al. [17, 18] presented a spectral two-and-a-half dimensional (2.5D) formulation of two elements to study guided waves in coupled problems involving thin-walled structures and30 fluid-acoustic enclosures. They proposed a plate element based on a Reissner-Mindlin/Kirchhoff-Love mixed formulation to represent thin-walled structures. These authors found accurate results using an appropriate combination of element size and order of the approximation functions. Regarding the development of the spectral boundary element method (SBEM), Calaon et al. [9] presented an approach to solve three dimensional (3D) elasticity problems using element shape functions generated35 from radial basis functions. They found that the exponential convergence of the solution was achievable. Holm et al. [10] extended the hp−version of the BEM to the Helmholtz equation to study the scattering of time harmonic acoustic waves at thin screens. Numerical experiments showed an exponential convergence rate and an arbitrarily high algebraic for the hp solution. Hwang [11] used a spectral approach to solve the boundary integral equation to study scattering and radiation problems. The spectral scheme presented40 orthogonal functions as basis functions to approximate the acoustic variables. Accurate results were found 2 with a reduction of computational effort with respect conventional BEM. Zou and Aliabadi [12] presented a boundary spectral element method for modelling high-frequency wave propagation in three dimensions. The authors used three types of high-order spectral elements for boundary discretisation. The computed results were validated against a well-stablished FEM. Milind and Mitra [15] studied time harmonic problems in45 acoustics using a Fourier-Hankel SBEM formulation. The authors compared the solution given by the SBEM with those obtained using the FEM, and other standard BEM techniques. Furthermore, Cerrato et al. [16] proposed a spectral coupled formulation of the BEM and the FEM to study linear water-wave propagation. They used the SFEM to study a closed region with non-constant bathymetry, and derived a SBEM from the mild-slope equation to study water-wave propagation over an infinite region with unidirectional variable50 bathymetry. They achieved very accurate results in comparison with the standard discrete solutions to this problem. In this work, we propose a two-and-a-half dimensional coupled spectral formulation based on the FEM and the BEM to study wave propagation in fluid and solid waveguides. The proposed method is useful for problems whose geometry and material properties are both uniform in one direction, whereas the excitation55 exhibits a three dimensional distribution. The spectral approach allows the analysis of wave propagation at high wavenumbers, which reduces the pollution effect with low computational effort. The outline of this work is as follows. First, the fluid-solid interaction problem is formulated in two- and-a-half dimensions. The 2.5D formulation of the BEM and the FEM are extended to their spectral formulations (SBEM and SFEM, respectively). The SBEM is used for the representation of both acoustic60 enclosures and unbounded domains. Moreover, a solid finite spectral element is proposed to represent elastic waveguides with arbitrary cross-section. The proposed approach is verified with two benchmark problems: the wave propagation in an unbounded acoustic region [16] and the scattering of waves by an elastic inclusion [19]. An h−pconvergence analysis is carried out to assess the accuracy of the method. Finally, the proposed method is applied to study the scattered wave field in an unbounded acoustic domain due to an array of65 elastic fluid-filled cylinders. 2. Numerical model The spectral formulations of the BEM and FEM are developed to represent an elastic waveguide submerged in an unbounded acoustic region, as shown in Figure 1. The proposed methodology uses spectral solid elements to represent the elastic waveguide (Ωs), while the limiting interface (Γ) between the solid and the fluid (Ωf) is modelled with spectral boundary elements. Both methods are coupled by imposing the appropriate boundary conditions at the interface Γ: equilibrium of normal pressure and compatibility of normal displacement must be ensured at the fluid-solid interface, with null shear stresses. The governing 3 Ωs Ωf Γ x y z Figure 1: Domain decomposition equations and the boundary conditions of the coupled fluid-solid system are:                ∇σ+b=ω2ρsusin Ωs,with BC ∇2pf+κ2 fpf= 0 in Ωf,with BC uT sn=ufin Γ, (σn)Tn+pf= 0 in Γ, (1) where the variables in the solid equation are the stress tensor σ, the body force vector b, the displacement vector us, the solid density ρsand the angular frequency ω. The Helmholtz equation defines the pressure field pffor a fluid wavenumber κf=ω/cf, where cfis the sound wave propagation velocity in the fluid70 subdomain. The coupling conditions include the fluid normal displacement to the boundary ufand the outward solid normal at the interface Γ, n. The 2.5D formulation is handled by defining a characteristic field for the elastic and the acoustic waveguides at point x=x(x, y, z) as [20, 21]: b f(x, κz, ω) = e f(e x, κz, ω)e−ικzz,(2) where e f(e x, κz, ω) is the frequency-wavenumber representation of an unknown variable (e.g. displacement, force or sound pressure), e x=e x(x, y) and ι=√−1. The 3D solution is computed as the superposition of 2D solutions with different wavenumber κz[22]: f(x, ω) = Z+∞ −∞ e f(e x, κz, ω)e−ικzzdκz.(3) The integral in Equation (3) presents poles and singularities, and the integration over the longitudinal wavenumber becomes mathematically and numerically complicated. One simple way around these difficulties is to replace the single-source problem, whose solution is expressed by Equation (3), with a multiple-source problem where sources are periodically distributed along the zaxis [22]. Then, Equation (3) is replaced by: 4 f(x, ω) = Z+∞ −∞ e f(e x, κz, ω)e−ικzz ∞ X m=−∞ eικmLdκz,(4) where Lis the periodicity source interval. Lneeds to be large enough to avoid the response being contaminated by the virtual source. Then, Equation (4) reduces to: f(x, ω) = 2π L ∞ X n=−∞ e f(e x, κzn, ω)e−ικznzdκz.(5) with κzn =2π Ln, which converges and can be approximated by a finite sum of Nterms as: f(x, ω) = 2π L N X n=−Ne f(e x, κzn, ω)e−ικznzdκz.(6) 2.1. Spectral boundary element method (SBEM) Next, the SBEM is derived from the standard collocation method, which is fully covered in the literature [23]. The main differences between the two methods concern the field variable approximation and the75 proposed element integration scheme. In both cases, the boundary integral representation is derived from the Helmholtz equation [23], considering a homogeneous fluid domain denoted by its volume Ωfand its boundary Γ. Although the proposed method is used to obtain the three-dimensional radiated field by an elastic waveguide with an arbitrary cross section, the problem is computed as the superposition of 2D problems with different wavenumber κz. Therefore, the boundary Γ is reduced to its cross section Σ at the80 coordinate plane z= 0. Thus, the integral representation of the fluid pressure for a point e xilocated at the boundary Γ in the frequency-wavenumber domain is: ci(e xi)epi(e xi, κz, ω) = −ZΣ ιρω2eui(e xi, κz, ω)e Ψ(e x, κz, ω;e xi) + epi(e xi, κz, ω)∂e Ψ(e x, κz, ω;e xi) ∂n!dΣ,(7) where ep(e xi, κz, ω) and eu(e xi, κz, ω) are the sound pressure and the particle normal displacement at the cross section Σ of the boundary Γ, respectively. e Ψ(e x, κz, ω;e xi) is the solution to the Helmholtz equation at point e xdue to a point source located at e xi[24]: e Ψ(e x, κz, ω;e xi) = −ι 4H(2) 0(kαr),(8) where r=ke x−e xikis the distance to the source, kα=qk2 f−κ2 zis the effective wavenumber and H(2) 0is the zero order Hankel function of second kind. The integral-free term ci(e xi) only depends on the boundary geometry at the collocation point e xi.85 The solution to this integral is found by discretisation of the boundary into Nelement with Σ = SN j=1 Σj: ci(e xi)epi(e xi, κz, ω) = − N X j=1 hehj i(e x, κz, ω;e xi)epj(e xj, κz, ω) + egj i(e x, κz, ω;e xi)euj(e xj, κz, ω)i,(9) 5 ξ -1 -0.5 0 0.5 1 -0.2 0 0.2 0.4 0.6 0.8 1φ1 φ2 φ3 φ4 φ5 φ6 φ7 Figure 2: Element shape functions φkfor an element with order p= 6. where egj iand ehj iare the boundary element matrices resulting from integration in element Σjand collocation at nodal point e xi. The field variables are approximated within the element using a set of interpolation shape functions φ(ξ) as: epj(e xj, κz, ω) = p+1 X k=1 φkep(e xk, κz, ω) = φepe, euj(e xj, κz, ω) = p+1 X k=1 φkeu(e xk, κz, ω) = φeue. (10) The proposed methodology uses Lagrange polynomials of order pas interpolation shape functions given by:90 φk=Y j6=i ξ−ξi ξj−ξi ,(11) where nodal coordinates ξiare found at the Lobatto-Gauss-Legendre (LGL) points. The use of a family of orthogonal polynomial (φi(ξj) = δij), such as Lagrange polynomials, as interpolation shape functions, mitigates the appearance of the Runge phenomenon, which can result in a highly ill-conditioned problem for high-order elements [25]. As an example, Figure 2 shows the shape functions for the element with order p= 6 used in this work.95 Once the sound pressure and particle normal displacement are approximated within the element, the former integrals in Equation (9) become: egj i(e x, κz, ω;e xi) = ZΣj φe Ψ (e x, κz, ω;e xi)dΣ,(12) ehj i(e x, κz, ω;e xi) = ZΣj φ∂e Ψ (e x, κz, ω;e xi) ∂n dΣ.(13) The spatial integrations in Equations (12) and (13) are numerically evaluated using a standard Gauss- Legendre quadrature with (p+ 1) integration points whenever the collocation point does not belong to 6 the integration element. However, the integrals become singular when the collocation point falls on the integration element. Then, the element integration is numerically evaluated by element subdivision, isolating the collocation point at the integration endpoints: ZΣje f(e x, κz, ω;e xi)dΣ = Z1 −1e f(e x, κz, ω;e xi)|Jj(ξ)|dξ =Zξp −1e f(e x, κz, ω;e xi)|Jj(ξ)|dξ +Z1 ξpe f(e x, κz, ω;e xi)|Jj(ξ)|dξ, (14) with ξpbeing the natural coordinate of the collocation point e xi,e f(e x, κz, ω;e xi) the kernel of integrals in Equations (12) and (13), and |Jj(ξ)|the Jacobian of the transformation to the natural coordinate system for element Σj. A domain transformation is used for the treatment of the singular integrals as follows [26]:    z2 l=ξp−ξ , −1≤ξ≤ξp z2 r=ξ−ξp, ξp≤ξ≤1 (15) Then, integrals in Equation (14) can be rewritten as: ZΣe f(e x, kz, ω;e xi)dΣ = Z0 √1+ξpe f(e x, kz, ω;e xi)J(ξ)dξ dzl dzl+Z√1−ξp 0e f(e x, kz, ω;e xi)J(ξ)dξ dzr dzr(16) where dξ dzland dξ dzrare respectively the Jacobian of the transformations obtained from Equation (15). The application of this polynomial transformation results in a significant improvement of convergence rates [26]. A linear transformation is used in order to transform the integration interval to η∈[−1,1]:    zl=√1+ξp 2(1 + ηl),−1≤ξ≤ξp zr=√1−ξp 2(1 + ηr), ξp≤ξ≤1 (17) and Equation (16) becomes: ZΣe f(e x, kz, ω;e xi)dΣ = Z1 −1e f(e x, kz, ω;e xi)J(ξ)dξ dzl dzl dηl dηl+Z1 −1e f(e x, kz, ω;e xi)J(ξ)dξ dzr dzr dηr dηr(18) where dzl dη and dzr dη are the Jacobian of the linear transformations appliad by Equation (17). In this work, the integrals in Equation (18) are numerically evaluated using a Gauss-Legendre quadrature with 10(p+ 1) integration points. The integration scheme allowed accurate results as can be seen in the following section.100 Finally, the boundary element matrices are assembled into a single system of equations relating sound pressure and particle normal displacement at nodal points: e H(e x, κz, ω)ep(e x, κz, ω) = e G(e x, κz, ω)eu(e x, κz, ω),(19) where ep(ex, κz, ω) and eu(e x, κz, ω) collect the nodal pressure and the normal displacement to the boundary, and e H(e x, κz, ω) and e G(e x, κz, ω) are the fully populated non-symmetrical boundary element system matrices. 7 Equation (19) can be modified in order to consider an incident pressure field to study scattering problems as: e H(e x, kz, ω)e p(e x, kz, ω) = e G(e x, kz, ω)e un(e x, kz, ω) + e pI(e x, kz, ω),(20) where e pIis an arbitrary incident pressure field. Once the boundary variables are solved, the radiated wave105 field at any receiver point in the fluid domain Γfis computed through the Somigliana identity. Then, the problem solution is defined as the superposition of the radiated wave field and the incident pressure field. The assembled boundary element matrices require special treatment if the boundary has normal discontinuities. Since the sound pressure is unique at any point on the boundary, the element assembling into the matrix e H(e x, κz, ω;e xi) is quite straightforward. However, the normal displacement at a discontinuity is not110 unique. This possibility is taken into account considering different values according to the element normals at this location. Nevertheless, the compatibility of normal displacements is not achieved when assembling matrix e G(e x, κz, ω;e xi). Then, the assembled matrix has one extra column for each discontinuity point at the boundary and the normal displacement vector also has an extra component. 2.2. Spectral finite element method (SFEM)115 The spectral finite element formulation is based on the virtual work principle in the frequency-wavenumber domain, which states [27]: −ω2ZΩs δb uT(x, κz, ω)ρsb u(x, κz, ω)dΩ + ZΩs δbεT(x, κz, ω)b σ(x, κz, ω)dΩ =ZΩs δb uT(x, κz, ω)ρsb b(x, κz, ω)dΩ + ZΓ δb uT(x, κz, ω)b q(x, κz, ω)dΓ, (21) where b u(x, κz, ω) is the displacement vector, bε(x, κz, ω) and b σ(x, κz, ω) are respectively the strain and stress tensors, ρsb b(x, κz, ω) is the body force in the domain Ωs,ρsis the solid density, and b q(x, κz, ω) is the surface force vector. A variable preceded by δdenotes a compatible variation of the displacement or the strain field. Voigh notation is used to write both the symmetrical stress tensor b σand the strain tensor bε, b σ={bσxx,bσyy,bσzz,bσxy,bσyz,bσzx}Tand bε={bεxx,bεyy,bεzz,bγxy,bγyz,bγzx}T, respectively.120 The stress tensor is related to the strain tensor through the constitutive relation: b σ=Cbε,(22) where, in the case of a linear isotropic material, the constitutive matrix Cdepends on the Young’s modulus 8 Eand the Poisson ratio ν: C=E (1 + ν)(1 −2ν)               1−ν ν ν 000 ν1−ν ν 000 ν ν 1−ν000 0001−2ν 20 0 0 0 0 0 1−2ν 20 0 0 0 0 0 1−2ν 2               .(23) The strain tensor can be derived from the displacement vector b uas: bε=Lb u,(24) where: L=               ∂ ∂x 0 0 0∂ ∂y 0 0 0 0 ∂ ∂y ∂ ∂x 0 0∂ ∂z ∂ ∂y ∂ ∂z 0∂ ∂x               .(25) The dependence of the strain tensor on the longitudinal coordinate zcan be removed considering the displacement definition in the frequency-wavenumber domain according to Equation (2): b u(x, κz, ω) = e u(e x, κz, ω)e−ικzz.(26) Then, Equation (24) can be rewritten as: bε=L1b u+L2 ∂b u ∂z =L1b u−ικzL2b u,(27) where the differential operators are: L1=               ∂ ∂x 0 0 0∂ ∂y 0 0 0 0 ∂ ∂y ∂ ∂x 0 0 0 ∂ ∂y 0 0 ∂ ∂x               ,L2=               000 000 001 000 010 100               .(28) The 2.5D formulation is derived assuming, as in the previous section, that solid material properties and geometry are both homogeneous in the longitudinal direction. Therefore, Equation (21) is further elaborated 9 Ωs Ωf Γ R e x0 epI x y z Figure 6: Elastic inclusion in an unbounded acoustic region 3.2. Wave scattering by a cylindrical elastic waveguide The proposed coupled spectral boundary and finite element method was verified from a benchmark problem involving a cylindrical elastic waveguide in a homogeneous unbounded fluid medium. The inclusion had a radius r= 0.6 m. The solid properties were compression wave propagation velocity cp= 2630 m/s, shear wave propagation velocity cs= 1416 m/s, and density ρs= 2250 kg/m3, while the fluid properties200 were sound wave propagation velocity cf= 1500 m/s and density ρf= 1000 kg/m3. The proposed SFEM represented the solid inclusion, while the SBEM was used to model the fluid domain. The elastic inclusion was subjected to a dilatational point source located in the fluid at e x0= (−3.5,0) from the cavity centre (Figure 6). The incident wave field epI(e x, κz, ω) was defined from the fundamental solution e Ψ(e x, κz, ω) [24] as: epI(e x, κz, ω)=2e Ψ(e x, κz, ω).(49) The problem solution was computed for a fluid wavenumber κf= 5.236 rad/m. A constant value of κz= 1.0 rad/m was assumed for comparison purposes. The analytical solution to this problem can be found in Reference [19]. The numerical results were compared with the reference solution and the scaled L2-error205 2was used to asses the accuracy of the proposed method. Figure 7 represents the real part of the analytical solution of the fluid pressure around the inclusion, and the horizontal and longitudinal solid displacements. Maximum amplitudes of the pressure are found near the source, because of the superposition of the incident and the scattered pressure fields, and behind the inclusion, due to the propagation of the incident wave through it. Both the pressure and the displacements210 present symmetry with respect the horizontal axis. Three discretisations were investigated. The characteristic element sizes were 1/h ={2.12, 4.24, 8.48}m−1, resulting in κfh={2.47,1.23,0.62}. The element order pwas set to obtain nodal densities per wavelength dλ={5,7.5,10,15,20}, which was limited to dλ= 20 for the finest discretisation because of the element size. Figure 8 shows a coarse, a medium and a fine problem discretisation that ensure dλ= 20.215 Table 1 summarizes the computed scaled L2-error 2and the CPU time for the different h−pdiscretisations used in this problem. The computed results presented a good agreement with the reference solution, 16 (a) (b) (c) Figure 7: Real part of the analytical solution for (a) total fluid pressure, (b) horizontal solid displacement and (c) longitudinal solid displacement. kf= 5.236 rad/m and kz= 1.0 rad/m as the error varied from 2≈1×10−3to 2≈1×10−5. However, the coarsest mesh did not accurately represent the problem solution for low element orders, since the nodal density per wavelength was not high enough. The smallest error was found for the medium size discretisation with an acceptable computational220 time, while the finest mesh had a higher error in a lower time. On the other hand, the coarsest mesh needed a higher computational effort to reach a similar level of error because of the higher element order approximation. Furthermore, the coarsest mesh presented an error decay up to an element order of p= 6, where the convergence was truncated. This behaviour comes from the integration of singularities in Equations (12) and (13). Then, the accuracy of the coupled method is limited by the accuracy of the SBEM.225 As previously commented, the CPU time in the SBEM depends on the element order and problem discretisation. Similarly, in the case of the SFEM, the necessary number of operations to do the element integrations was nelem(p+ 1)2for a given mesh with nelem elements. The CPU time of the SFEM was proportional to the total number of operations depending on the element order by C0(p). Then, the CPU 17 x [m] -0.6 -0.3 0 0.3 0.6 y [m] -0.6 -0.3 0 0.3 0.6 (a) x [m] -0.6 -0.3 0 0.3 0.6 y [m] -0.6 -0.3 0 0.3 0.6 (b) x [m] -0.6 -0.3 0 0.3 0.6 y [m] -0.6 -0.3 0 0.3 0.6 (c) Figure 8: Problem discretisation using (a) coarse, (b) medium, and (c) fine meshes with nodal density dλ= 20. SBEM nodes are represented by circles and SFEM nodes by points time can be described as C0(p)nelem(p+1)2. Moreover, the computational cost of the SFEM was much more230 lower than the SBEM while, in comparison, the CPU time of coupling and solving the system of equations was practically negligible. Again, a hp-refinement is more efficient in terms of error and computational effort than a single h-refinement or p-refinement. 4. Numerical example Finally, the capabilities of the proposed method are explored with a numerical example. The problem235 solved herein was the study of the scattered wave field by an array of elastic cylinders submerged in an inviscid fluid medium (Figure 9). The elastic cylinders had an external radius r= 0.4 m and a thickness t= 0.05 m. The centres of the scatterers were equally spaced s= 1.2 m. The material was assumed to be PVC, with compression wave propagation velocity cp= 2143 m/s, shear wave propagation velocity cs= 875 m/s, and density ρ= 1400 kg/m3. The outer fluid was water, with sound wave propagation240 velocity cf= 1500 m/s and density ρf= 1000 kg/m3. Three configurations were considered based on the interior of the scatterers: i) water-filled; ii) air-filled, with sound propagation velocity 340 m/s and density 1.225 kg/m3; and iii) state of vacuum. The array of scatterers was subjected to a dilatational point source located in the fluid at position e x0= (−4.0,0) (Figure 9). The incident wave field epIwas defined as in the previous section. The problem245 response of the fluid-structure system was studied for a frequency f= 6400 Hz. The scatterers were modelled with spectral finite elements while outer and inner fluid subdomains were represented with the SBEM. The discretisation was chosen to enable a nodal density per wavelength of dλ= 12 in the outer fluid, while the scatterers and the inner fluid meshes were defined to match the outer fluid mesh. The outer fluid was modelled with 34 elements, so the cylinders and the inner fluid were modelled250 18 Table 1: Summary of the number of elements, nelem, number of nodes nnode for both SBEM and SFEM, scaled L2-error 2of the total fluid pressure, horizontal and longitudinal solid displacement and CPU time, computed for the different hp discretisations. SBEM SFEM Scaled L2error 2CPU time [s] 1/h[m−1]nelem nnode nelem nnode p dλp 2ux 2uz 2 16 57 2 5 2.1362 ×10−12.2429 ×10−15.7590 ×10−20.526 24 121 3 7.5 1.2179 ×10−21.5760 ×10−28.0715 ×10−31.051 2.12 8 32 12 209 4 10 1.7836 ×10−32.0796 ×10−31.0771 ×10−31.713 48 457 6 15 5.7194 ×10−44.9001 ×10−41.0273 ×10−34.379 64 801 8 20 5.0114 ×10−43.9949 ×10−41.0014 ×10−411.36 32 209 2 10 3.9259 ×10−35.4951 ×10−33.9152 ×10−30.899 4.24 16 48 48 457 3 15 1.2872 ×10−43.0848 ×10−42.3735 ×10−42.018 64 801 4 20 6.0030 ×10−52.7535 ×10−57.5524 ×10−53.976 8.48 32 64 192 801 2 20 1.6663 ×10−44.2036 ×10−43.2369 ×10−42.449 Ωs Ωf Ωf Γ Γ R s s s/2 e x0 epIx y z Figure 9: Configuration of scatterers. 19 x[m] -0.4 -0.2 0 0.2 0.4 y[m] -0.4 -0.2 0 0.2 0.4 Figure 10: Scatter discretisation with dλ= 12 and κfh= 2 again with 34 elements. The element order was set to p= 4 in all cases. This discretisation resulted in a nodal density dλ= 12 with and κfh= 2 in the outer fluid. Figure 10 shows the discretisation of one scatter. The longitudinal wavenumber κzwas not taken into account when computing the nodal density per wavelength. The 3D solution was computed as the superposition of 2D problems with different longitudinal wavenumbers by means of Equation (6). The 3D solution can be computed in a set of regularly spaced points along the longitudinal direction where ∆z= 2π/max {κz}. The maximum wavenumber max {κz}= 160 rad/m was chosen to compute the solution enabling six points in a longitudinal wavelength, ∆z= 2π/max {κz}= 0.0393 m. The maximum wavenumber max {κz}should be large enough to account for the slowest waves, namely guided waves, given by: vf=2πf max {κz}.(50) In this case, −160 rad/m≤κz≤160 rad/m, so it is possible to account for waves with velocity higher than255 250 m/s. The wavenumber sampling ∆κz= 0.5178 rad/m allowed to compute the solution for a maximum distance L=π/∆κz= 6.0672 m from the source. First, the behaviour of the cross section is analysed. Figure 11 shows the total pressure field in the outer fluid and the scattered pressure field inside the cylinders at a vertical plane z= 0 for the three studied cases. The total pressure fields are quite similar in the three cases. The maximum pressure was found at the260 source point and the scatterer configurations produced a shadowed zone, with a lower pressure amplitude. The vibration of the internal surface of the cylinders caused a scattered wave field inside them in the fluidfilled cases. When the cylinders were filled with water (Figure 9(a)), the scattered pressure field shows a distribution in accordance with the incident field, similar wavelength, amplitude, phase and direction, as the fluid properties were the same. The fluid-filled case (Figure 9(b)) exhibits a much more complex pattern,265 20 (a) (b) (c) Figure 11: Real part of the pressure field at cross-section of the scatterers system at the vertical plane z= 0 for the (a) water-filled cylinders, (b) air-filled cylinders and (c) vacuum cylinders. Sound pressure in water medium was normalised to 0.15 Pa, and to 5 ×10−5Pa in air medium. with lower amplitude and wavelength and any concordance with the incident field in direction or phase due 21 Figure 12: Real part of the pressure field (normalised to 0.15 Pa) in an unbounded fluid by a system of water-filled elastic scatterers. to the difference in the properties of inner and outer fluids. The vacuum and air-filled systems presented the best attenuation, as the water-filled system allowed the waves to propagate through the scatterers. Figure 12 shows the 3D representation of the sound pressure at the surface of the scatterers and the total pressure field in the fluid domain for the water-filled example. The scattered pressure field at a vertical plane270 z= 6.0672 m inside the cylinders is also shown. Only half of the model is represented for clarity’s sake, because of problem symmetry. The pressure field had a spherical distribution centred at the source, which was progressively decreasing as the distance to the source increased. Figure 13 shows the deformed shape of the solid scatterers external surface and displays the normal displacement to the boundary in the colour scale. The displacement field also had a spherical distribution centered at the source, which decreased as275 the distance from the source increased. Maximum amplitudes of both the pressure and displacement were found at the plane z= 0, where the source was located and the incident field was higher. This example has shown that the proposed method can be used to study the wave propagation in a fluid-filled elastic waveguides immersed in an unbounded acoustic medium. The proposed method was able to represent different fluid properties for the unbounded fluid and the acoustic enclosures.280 22 Figure 13: Real part of normal displacement at external surface (normalised to 1.2×10−12 m) in a system of water-filled elastic scatterers, represented over the deformed shape. 23 5. Conclusions This work has proposed a spectral element formulation based on the BEM and the FEM to study fluid and solid wave propagation. The presented methodology looks at 3D problems whose materials and geometric properties remain homogeneous in one direction. Solid subdomains were modelled with the SFEM, whereas the unbounded fluid media and acoustic enclosures were represented by the SBEM. The coupling285 of both methods was carried out by imposing the appropriate boundary conditions at the interface to study fluid-structure wave propagation. The model was verified with two benchmark problems with analytical solutions. The SBEM was verified with a problem involving wave propagation over an unbounded acoustic domain. The coupled SBEM-SFEM was verified with a problem concerning a cylindrical elastic waveguide in an unbounded fluid medium.290 Numerical results show good agreement with the analytical solution. An h−panalysis shows that an optimal solution in terms of accuracy and computational effort can be obtained through an hp-refinement. Discretisations with a nodal density per wavelength of ten are enough to achieve reasonable accuracy. The accuracy improves by increasing the nodal density up to twenty. To conclude, we found that medium size discretisations gave as accurate results as coarse and fine meshes for the same nodal density per wavelenght,295 with lower computational cost. Finally, the coupled methodology was used to study the wave field radiated by a three dimensional fluid-filled elastic scatterers system submerged in an unbounded acoustic medium. The sound pressure field and the scatterer’s surface displacements were studied. The proposed method made it possible to represent different fluids for the unbounded medium and the acoustic enclosures.300 Acknowledgements The research work presented herein was supported by the Spanish Ministry of the Economy and Competitiveness (Ministerio de Econom´ıa y Competitividad) through research project BIA2016-75042-C2-1-R. It was also supported by project POCI-01-0247-FEDER-017759 (SmartCore – Desenvolvimento de pain´eis multifuncionais de elevada performance), funded by the Operational Program for Competitiveness and In-305 ternationalization (POCI) of Portugal 2020, with the support of the European Regional Development Fund (FEDER). Financial support is gratefully acknowledged. The authors also wish to acknowledge the support provided by the Andalusian Scientific Computing Centre (CICA). References [1] I. Babuˇska, F. Ihlenburg, E. T. Paik, S. A. Sauter, A Generalized Finite Element Method for solving the Helmholtz310 equation in two dimensions with minimal pollution, Computer Methods in Applied Mechanics and Engineering 128 (3) (1995) 325 – 359. 24 [2] F. Ihlenburg, I. Babuˇska, S. Sauter, Reliability of finite element methods for the numerical computation of waves, Advances in Engineering Software 28 (7) (1997) 417 – 424. [3] A. A. Oberai, P. M. Pinsky, A numerical comparison of finite element methods for the Helmholtz equation, Journal of315 Computational Acoustics 08 (01) (2000) 211–221. [4] C. Willberg, S. Duczek, J.M. Vivar Perez, D. Schmicker, U. Gabbert, Comparison of different higher order finite element schemes for the simulation of Lamb waves, Computer Methods in Applied Mechanics and Engineering 241-244 (2012) 246 – 261. [5] A. T. Patera, A spectral element method for fluid dynamics: Laminar flow in a channel expansion, Journal of Computa-320 tional Physics 54 (3) (1984) 468 – 488. [6] W. Ostachowicz, P. Kudela, M. Krawczuk, A. Zak, Guided Waves in Structures for SHM: The Time - Domain Spectral Element Method, Wiley, 2012. [7] G. Degrande, G. De Roeck, A spectral element method for two-dimensional wave propagation in horizontally layered saturated porous media, Computers & Structures 44 (4) (1992) 717 – 728, Special Issue: Computational Structures325 Technology. [8] P. Z. Bar-Yoseph, D. Fisher, O. Gottlieb, Spectral element methods for nonlinear spatio-temporal dynamics of an Euler- Bernoulli beam, Computational Mechanics 19 (1) (1996) 136 – 151. [9] A. Calaon, R. Adey, J. Baynham, A spectral boundary element method, in: A. Kassab, C. A. Brebbia, E. Divo, D. Poljak (Eds.), Boundary Elements XXVII, Vol. 39, WIT Press, 2005, p. 65.330 [10] H. Holm, M. Maischak, E.P. Stephan, The hp-Version of the Boundary Element Method for Helmholtz Screen Problems, Computing, 57 (2) (1996) 105–134. [11] W.S. Hwang, Boundary spectral method for acoustic scattering and radiation problems, Journal of the Acoustical Society of America, 102(1) (1997) 96–101. [12] F. Zou, M.H. Aliabadi, On modelling three-dimensional elastodynamic wave propagation with boundary spectral element335 method, European Journal of Computational Mechanics, 27(3) (2018) 204–228. [13] P. Kudela, M. Krawczuk, W. Ostachowicz, Wave propagation modelling in 1D structures using spectral finite elements, Journal of Sound and Vibration 300 (1) (2007) 88 – 100. [14] C. Zhu, G. Qin, J. Zhang, Implicit Chebyshev spectral element method for acoustics wave equations, Finite Elements in Analysis and Design 47 (2) (2011) 184 – 194.340 [15] T. R. Milind, M. Mitra, A Study on Fourier Spectral and Fourier-Hankel Spectral Field Variable Approximations in BEM for 2D Interior Helmholtz Problems, Procedia Engineering 144 (2016) 520 – 527, International Conference on Vibration Problems 2015. [16] A. Cerrato, L. Rodr´ıguez-Tembleque, J. A. Gonz´alez, M. H. Ferri Aliabadi, A coupled finite and boundary spectral element method for linear water-wave propagation problems, Applied Mathematical Modelling 48 (2017) 1 – 20.345 [17] A. Romero, F.J. Cruz-Mu˜noz, A. Tadeu, P. Galv´ın, A novel 2.5D spectral approach for studying thin-walled waveguides with fluid-acoustic interaction, Computers and Structures 204 (2018) 1 – 19. [18] F.J. Cruz-Mu˜noz, A. Romero, A. Tadeu, and P. Galv´ın, A 2.5D spectral approach to represent acoustic and elastic waveguides interaction on thin slab structures, Procedia Engineering 199 (2017) 1374–1379. [19] A. Pereira, A. Tadeu, J. Ant´onio, Influence of the cross-section geometry of a cylindrical solid submerged in an acoustic350 medium on wave propagation, Wave Motion 36 (1) (2002) 23–39. [20] L. Gavri´c, Computation of propagative waves in free rail using a finite element technique, Journal of Sound and Vibration 185 (3) (1995) 531 – 543. [21] L. Gavri´c, Finite element computation of dispersion properties of thin-walled waveguides, Journal of Sound and Vibration 173 (1) (1994) 113–124.355 25