A 2.5D BEM-FEM using a spectral approach to study scattered waves in fluid–solid interaction problems

A 2.5D BEM-FEM USING A SPECTRAL APPROACH TO STUDY SCATTERED WAVES IN FLUID–SOLID INTERACTION PROBLEMS FRANCISCO JAVIER CRUZ-MU ˜ NOZ1, ANTONIO ROMERO1, PEDRO GALV´ IN1& ANT ´ ONIO TADEU2,3 1Departamento de Mec´ anica de Medios Continuos y Teor´ ıa de Estructuras, Escuela T´ ecnica Superior de Ingenier´ ıa, Universidad de Sevilla, Spain 2ITeCons – Institute for Research and Technological Development in Construction, Energy, Environment and Sustainability, Portugal 3ADAI – LAETA, Department of Civil Engineering, University of 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 a benchmark problem regarding 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 frequency range, the standard formulation using linear elements accurately represents the fluid and solid scattering waves. However, at high frequencies, the solution deteriorates due to so-called pollution effects [1], [2]. Higher element 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 because of the necessary computational effort to represent small wavelengths [4]. Many works suggest using high-order 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. 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 excitation exhibits a three dimensional distribution. www.witpress.com, ISSN 1743-3533 (on-line) WIT Transactions on Engineering Sciences, Vol 126, © 2019 WIT Press doi:10.2495/BE420101 Boundary Elements and other Mesh Reduction Methods XLII 111 Ωs Ωf Γ x y z Figure 1: Domain decomposition. 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. The2.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 acoustic 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 [7] and the scattering of waves by an elastic inclusion [8]. 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 of 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 Fig. 1. The proposed methodology uses spectral solid elements to represent the elastic waveguide (Ωs), while the fluid domain (Ω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 fluidsolid interface, with null shear stresses. 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 [9], [10]: b f(x, κz, ω) = e f(ex, κz, ω)e−ικzz,(1) where e f(ex, κz, ω)is the frequency-wavenumber representation of an unknown variable (e.g. displacement, force or sound pressure), ex=ex(x, y)and ι=√−1. The 3D solution is computed as the superposition of 2D solutions with different wavenumber κz: f(x, ω) = Z+∞ −∞ e f(ex, κz, ω)e−ικzzdκz.(2) 2.1 Spectral boundary element method (SBEM) Next, the SBEM is derived from the standard collocation method, which is fully covered in the literature [11]. The main differences between the two methods concern the field Figure 1: Domain decomposition. 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 acoustic 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 [7] and the scattering of waves by an elastic inclusion [8]. 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 of 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 Fig. 1. The proposed methodology uses spectral solid elements to represent the elastic waveguide (Ωs), while the fluid domain (Ω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 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 [9], [10]: b f(x, κz, ω) = e f(ex, κz, ω)e−ικzz,(1) where e f(ex, κz, ω)is the frequency-wavenumber representation of an unknown variable (e.g. displacement, force or sound pressure), ex=ex(x, y)and ι=√−1. The 3D solution is computed as the superposition of 2D solutions with different wavenumber κz: f(x, ω) = Z+∞ −∞ e f(ex, κz, ω)e−ικzzdκz.(2) www.witpress.com, ISSN 1743-3533 (on-line) WIT Transactions on Engineering Sciences, Vol 126, © 2019 WIT Press 112 Boundary Elements and other Mesh Reduction Methods XLII 2.1 Spectral boundary element method (SBEM) Next, the SBEM is derived from the standard collocation method, which is fully covered in the literature [11]. The main differences between the two methods concern the field variable approximation and the proposed element integration scheme. In both cases, the boundary integral representation is derived from the Helmholtz equation [11], considering a homogeneous fluid domain denoted by its volume Ωfand 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 the coordinate plane z= 0. Thus, the integral representation of the fluid pressure for a point exi located at the boundary Γin the frequency-wavenumber domain is: ci(exi)epi(exi, κz, ω) =−ZΣ ιρω2eui(exi, κz, ω)e Ψ(ex, κz, ω;exi) + epi(exi, κz, ω)∂e Ψ(ex, κz, ω;exi) ∂n!dΣ, (3) where ep(exi, κz, ω)and eu(exi, κz, ω)are the sound pressure and the particle normal displacement at the cross section Σof the boundary Γ, respectively. e Ψ(ex, κz, ω;exi)is the solution to the Helmholtz equation at point exdue to a point source located at exi[12]. The integral-free term ci(exi)only depends on the boundary geometry at the collocation point exi. The boundary is discretised into Nelements, with Σ = SN j=1 Σj, and the field variables are approximated within the element using a set of interpolation shape functions φ(ξ). The proposed methodology uses Lagrange polynomials of order pas interpolation shape functions given by: φk=Y j6=i ξ−ξi ξj−ξi ,(4) 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 [13]. Once the sound pressure and particle normal displacement are approximated within the element, the former integrals in eqn (3) become: egj i(ex, κz, ω;exi) = ZΣj φe Ψ (ex, κz, ω;exi)dΣ,(5) ehj i(ex, κz, ω;exi) = ZΣj φ∂e Ψ (ex, κz, ω;exi) ∂n dΣ.(6) The spatial integrations in eqns (5) and (6) are numerically evaluated using a standard Gauss- Legendre quadrature with (p+ 1) integration points whenever the collocation point does not belong to the integration element. However, the integrals become singular when the collocation point falls on the integration element. These integrals are numerically evaluated www.witpress.com, ISSN 1743-3533 (on-line) WIT Transactions on Engineering Sciences, Vol 126, © 2019 WIT Press Boundary Elements and other Mesh Reduction Methods XLII 113 by element subdivision, isolating the collocation point at the integration endpoints ZΣje f(ex, κz, ω;exi)dΣ = Z1 −1e f(ex, κz, ω;exi)|Jj(ξ)|dξ =Zξp −1e f(ex, κz, ω;exi)|Jj(ξ)|dξ +Z1 ξpe f(ex, κz, ω;exi)|Jj(ξ)|dξ, (7) with ξpbeing the natural coordinate of the collocation point exi,e f(ex, κz, ω;exi)the kernel of integrals in eqns (5) and (6), and |Jj(ξ)|the Jacobian of the transformation to the natural coordinate system for element Σj. These integrals are numerically evaluated using a smooth polynomial transformation as described in Reference [14]. Accurate results are given in the next section for a number of integration points 10(p+ 1). 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(ex, κz, ω)ep(ex, κz, ω) = e G(ex, κz, ω)eu(ex, κz, ω),(8) where ep(ex, κz, ω)and eu(ex, κz, ω)collect the nodal pressure and the normal displacement to the boundary, and e H(ex, κz, ω)and e G(ex, κz, ω)are the fully populated non-symmetrical boundary element system matrices. 2.2 Spectral finite element method (SFEM) The spectral finite element formulation is based on the virtual work principle in the frequencywavenumber domain, which states [15]: −ω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Γ,(9) 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 nodal force vector. A variable preceded by δdenotes a compatible variation of the displacement or the strain field. The stress tensor is related to the strain tensor through the constitutive relation b σ=Cbε, where, in the case of a linear isotropic material, the constitutive matrix Cdepends on the Young’s modulus Eand the Poisson ratio ν. The dependence of the strain tensor on the longitudinal coordinate zcan be removed considering the displacement definition in the frequency-wavenumber domain according to eqn (1): b u(x, κz, ω) = e u(ex, κz, ω)e−ικzz.(10) Then, the strain tensor can be derived from the displacement vector b uas: bε=L1b u+L2 ∂b u ∂z =L1b u−ικzL2b u,(11) where the differential operators are: www.witpress.com, ISSN 1743-3533 (on-line) WIT Transactions on Engineering Sciences, Vol 126, © 2019 WIT Press 114 Boundary Elements and other Mesh Reduction Methods XLII 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        .(12) The 2.5D formulation is derived assuming that solid material properties and geometry are both homogeneous in the longitudinal direction. Therefore, eqn (9) is further elaborated by eliminating the dependency on the coordinate zrewriting the volume integrals over the section Asof the domain integration Ωsby the coordinate plane z= 0: −ω2ZAs δe uTρse udA +ZAs δeεTe σdA =ZAs δe uTρse bdA +ZΣ δe uTe qdΣ,(13) where Σis defined as the boundary Γat z= 0. The solid subdomain is discretised into elements, where the approximated displacement vector e uis defined as: e u=ϕe ue,(14) where ϕare the two-dimensional shape functions obtained from the one-dimensional function φ(ξ)as: ϕk(ξ, η) = φi(ξ)φj(η), k = (i−1)(p+ 1) + j. (15) The proposed spectral finite element is defined from local nodal coordinates (ξ, η)∈ ([−1,1] ×[−1,1]) at the LGL points. Once the element approximation is defined, the strain tensor is expressed in terms of the nodal displacement as: eεe=L1ϕe ue−ικzL2ϕe ue=B1e ue−ικzB2e ue,(16) where B1=L1ϕand B2=L2ϕ. Introducing eqns (14) and (16) into eqn (13) and operating, the virtual work principle can be rewritten for each element as: −ω2ZAe s δe ueT ϕTρsϕe uedA +ZAe s δe ueT BT 1CB1e uedA −ικzZAe s δe ueT BT 1CB2+BT 2CB1e uedA −κ2 zZAe s δe ueT BT 2CB2e uedA =ZAe s δe ueT ϕTρse bedA +ZΣe δe ueT ϕTe qedΣ,(17) where Ae sis the cross-sectional area of the three dimensional element in the zdirection and Σeits boundary. The following expression is obtained taking into account that eqn (17) must be satisfied for any compatible displacement δe ue: −ω2Me+Ke 0−ικzKe 1−κ2 zKe 2e ue=efe,(18) www.witpress.com, ISSN 1743-3533 (on-line) WIT Transactions on Engineering Sciences, Vol 126, © 2019 WIT Press Boundary Elements and other Mesh Reduction Methods XLII 115 where Meis the element mass matrix, Ke 0,Ke 1and Ke 2are the stiffness element matrices and bfeis a vector that collects the forces applied to the element. These matrices are defined by: Me=ZAe s ϕTρsϕdA, (19) Ke 0=ZAe s BT 1CB1dA, (20) Ke 1=ZAe sBT 1CB2−BT 2CB1dA, (21) Ke 2=ZAe s BT 2CB2dA, (22) and the external load vector efeas: efe(ex, κz, ω) = ZAe s ϕTρse be(ex, κz, ω)dA +ZΣe ϕTe qe(ex, κz, ω)dΣ,(23) which contains the contributions of both body forces and surface loads. Then, eqn (18) can be written as: e Ke u=ef,(24) where e K= [−ω2M+K0−ικzK1−κ2 zK2]is the dynamic stiffness matrix obtained in the frequency-wavenumber domain from the assembled mass Mmatrix, and stiffness K0,K1and K2matrices. The mass Mematrix and stiffness Ke 0,Ke 1and Ke 2matrices are numerically integrated using a LGL quadrature of order p. Therefore, the nodal element coordinates coincide with the element integration points producing a diagonal mass Mematrix and stiffness Ke 2matrices. 2.3 SBEM-SFEM coupling Eqns (8) and (24) are coupled by imposing equilibrium of forces, compatibility of normal displacement and null shear stresses at the interface Γ. These equations are assembled into a single comprehensive system, together with the equilibrium and compatibility conditions. The load vector effat the fluid–solid interface is obtained by integrating the fluid pressure field epfover the boundary Γ: eff=−ZΓ ϕTnφepfdΓ = −RTepf,(25) where nis the outward normal vector along Γ, and Ris the coupling fluid–solid matrix which relates force at the solid subdomain and pressure at the boundary. Substituting eqn (25) into eqn (24) yields: e Kss e Ksf e Kfs e Kff e us e uf=efs −RTepf,(26) where subscript findicates degrees of freedom belonging to Γand sstands for the rest of solid degrees of freedom. Finally, the coupling of eqns (8) and (26) is carried out with the imposition of compatibility of normal displacement at the boundary Γ, and null shear-stresses. Both www.witpress.com, ISSN 1743-3533 (on-line) WIT Transactions on Engineering Sciences, Vol 126, © 2019 WIT Press 116 Boundary Elements and other Mesh Reduction Methods XLII Ωs Ωf Γ R e x0 epI x y z Figure 2: Elastic inclusion in an unbounded acoustic region. systems of equations are assembled into an overall system:  e Kss e Ksf 0 e Kfs e Kff RT 0−e GNTe H  e us e uf epf = efs eff 0 ,(27) where Nis a matrix containing the outward normal vector nat the boundary. Eqn (27) is solved for each frequency-wavenumber step to compute the coupled fluid-solid response. 3 MODEL VERIFICATION The proposed coupled spectral boundary and finite element method was verified on 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 properties were sound wave propagation velocity cf= 1500 m/sand 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 distance ex0= (−3.5,0) from the cavity centre (Fig. 2). The problem solution was computed for a fluid wavenumber κf= 5.236 rad/mover a grid of receivers using a density of six points per wavelength. A constant value of κz= 1.0 rad/mwas assumed for comparison purposes. The analytical solution to this problem can be found in Reference [8]. The numerical results were compared with the reference solution and the scaled L2-error 2was used to asses the accuracy of the proposed method. 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λ={10,15,20}, which was limited to dλ= 20 for the finest discretisation because of the element size. 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, as the error varied from 2≈1×10−3to 2≈1×10−5. The smallest error was found for the medium size discretisation with an acceptable computational time, while the finest mesh had a higher error in a similar time. On the other hand, the coarsest mesh needed a higher computational effort to reach an acceptable error because of the higher element order approximation. Therefore, a hp-refinement is more efficient in terms of error and computational effort than a single h-refinement or p-refinement. Figure 2: Elastic inclusion in an unbounded acoustic region. systems of equations are assembled into an overall system:  e Kss e Ksf 0 e Kfs e Kff RT 0−e GNTe H  e us e uf epf = efs eff 0 ,(27) where Nis a matrix containing the outward normal vector nat the boundary. Eqn (27) is solved for each frequency-wavenumber step to compute the coupled fluid–solid response. 3 MODEL VERIFICATION The proposed coupled spectral boundary and finite element method was verified on 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 properties were sound wave propagation velocity cf= 1500 m/sand 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 distance ex0= (−3.5,0) from the cavity centre (Fig. 2). The problem solution was computed for a fluid wavenumber κf= 5.236 rad/mover a grid of receivers using a density of six points per wavelength. A constant value of κz= 1.0 rad/mwas assumed for comparison purposes. The analytical solution to this problem can be found in Reference [8]. The numerical results were compared with the reference solution and the scaled L2-error 2was used to asses the accuracy of the proposed method. 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λ={10,15,20}, which was limited to dλ= 20 for the finest discretisation because of the element size. 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, as the error varied from 2≈1×10−3to 2≈1×10−5. The smallest error was found for the medium size discretisation with an acceptable computational www.witpress.com, ISSN 1743-3533 (on-line) WIT Transactions on Engineering Sciences, Vol 126, © 2019 WIT Press Boundary Elements and other Mesh Reduction Methods XLII 117 Table 1: Summary of scaled L2-error 2of the total fluid pressure, horizontal and longitudinal solid displacement and CPU time, computed for the different hp discretisations. 1/h[m−1]p dλp 2ux 2uz 2CPU time [s] 2 5 2.14 ×10−12.24 ×10−15.76 ×10−20.526 3 7.5 1.22 ×10−21.58 ×10−28.07 ×10−31.051 2.12 4 10 1.78 ×10−32.08 ×10−31.08 ×10−31.713 6 15 5.71 ×10−44.90 ×10−41.03 ×10−34.379 8 20 5.01 ×10−43.99 ×10−41.00 ×10−411.36 2 10 3.93 ×10−35.50 ×10−33.92 ×10−30.899 4.24 3 15 1.29 ×10−43.08 ×10−42.37 ×10−42.018 4 20 6.0×10−52.75 ×10−57.55 ×10−53.976 8.48 2 20 1.67 ×10−44.20 ×10−43.24 ×10−42.449 time, while the finest mesh had a higher error in a similar time. On the other hand, the coarsest mesh needed a higher computational effort to reach an acceptable error because of the higher element order approximation. Therefore, 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 problem solved herein was a study of the wave field scattered by an array of elastic cylinders submerged in an inviscid fluid medium (Fig. 3). 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 propagation velocity cf= 1500 m/sand 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/sand 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 ex0= (−4.0,0) (Fig. 3). The incident wave field epIwas defined as in the previous section. The problem 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 modelled 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. The 3D solution was computed as the superposition of 2D problems with different longitudinal wavenumbers by means of eqn (2). 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/mwas chosen to compute the solution enabling six points in a longitudinal wavelength, ∆z= 2π/max {κz}= 0.0393 m. The www.witpress.com, ISSN 1743-3533 (on-line) WIT Transactions on Engineering Sciences, Vol 126, © 2019 WIT Press 118 Boundary Elements and other Mesh Reduction Methods XLII Ωs Ωf Ωf Γ Γ R s s s/2 e x0 epIx y z Figure 3: Configuration of scatterers. maximum pressure was found at the 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 fluid-filled cases. When the cylinders were filled with water (Fig. 4(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 (Fig. 4(b)) exhibited a much more complex pattern, with lower amplitude and wavelength and any concordance with the incident field in direction or phase due 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. Fig. 5 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 plane 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. Fig. 6 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 as 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. Figure 3: Configuration of scatterers. wavenumber sampling ∆κz= 0.5178 rad/mmade it possible to compute the solution for a maximum distance L=π/∆κz= 6.0672 m from the source. First, the behaviour of a cross section is analysed. Fig. 4 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 the 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 fluid-filled cases. When the cylinders were filled with water (Fig. 4(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 (Fig. 4(b)) exhibited a much more complex pattern, with lower amplitude and wavelength and any concordance with the incident field in direction or phase due 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. Fig. 5 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 plane 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. Fig. 6 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 www.witpress.com, ISSN 1743-3533 (on-line) WIT Transactions on Engineering Sciences, Vol 126, © 2019 WIT Press Boundary Elements and other Mesh Reduction Methods XLII 119