A Hybrid Analytical-Numerical Model Based on the Method of Fundamental Solutions for the Analysis of Sound Scattering by Buried Shell Structures

Several numerical and analytical models have been used to study underwater acoustics problems. The most accurate and realistic models are usually based on the solution of the wave equation using a variety of methods. Here, a hybrid numerical-analytical model is proposed to address the problem of underwater sound scattering by an elastic shell structure, which is assumed to be circular and that is buried in a fluid seabed bellow a water waveguide. The interior of the shell is filled with a fluid that may have different properties from the host medium. The analysis is performed by coupling analytical solutions developed both for sound propagation in thewaveguide and in the vicinity of the circular hollow pipeline. The coupling between solutions is performed using the method of fundamental solutions. This strategy allows a compact description of the propagation medium while being very accurate and highly efficient from the computational point of view.


Introduction
The detection of buried objects in solid and fluid media has been an active research topic, making use of different approaches.Techniques based on wave propagation have received particular interest from researchers, leading to the development of a broad variety of analytical and numerical models to simulate this propagation and to an intense research on the interpretation of field results.In fact, the measurement of spatial and temporal variations, recorded at hydrophones or geophones, resulting from the generation of waves produced by dynamic sources, placed inside elastic or acoustic media, is frequently used to infer the presence of buried structures or the geological configuration of a specific site.Some early reference works on this topic are due to Claerbout 1 or Griffiths and King 2 , addressing the specific application of such methods to geophysics; in the case of underwater acoustics, reference works include the now classic book by Jensen et al. 3 , which describes a number of formulations that can be used in both shallow water and deep water scattering problems.
In this paper, the authors focus on the scattering of waves in underwater configurations, for which different methods have been used in the past, ranging from the analytical methods presented by Pao and Mow 4 for studying wave diffraction near cylindrical circular inclusions to purely numerical methods, such as finite difference e.g., Stephens 5 or Wang 6 and finite elements techniques e.g., Marfurt 7 or Zampolli et al. 8 , combined with transmitting boundaries.
Methods relying on the description of the relevant boundaries of the problem have also been developed and form a very interesting class for this type of applications.An important early work on acoustic scattering in the open ocean using the boundary element method is due to Dawson and Fawcett 9 , which analyses a waveguide with flat surfaces, except for a compact area of deformation, where the acoustic scattering takes place.A hybrid model which combines the standard boundary element method BEM in an inner region with varying bathimetry and an eigenfunction expansion in the outer region of constant depth was later proposed by Grilli et al. 10 .Works by Santiago and Wrobel 11 described the implementation of the subregion technique in boundary element formulation for the analysis of two-dimensional acoustic wave propagation in a shallow water region with irregular seabed topography, allowing for the analysis of more general underwater systems.In their approach, the bottom and surface boundaries of the regions are modeled using Neumann and Dirichlet conditions, allowing for the use of Green's functions that satisfy either the free surface boundary condition or both the boundary conditions on the free surface and rigid bottom.
The use of specific Green's functions that account for part of the boundaries of the analysis domain has been an important strategy when dealing with boundary element methods, since they may allow for smaller discretization schemes, leading to lower computational effort, and therefore many researchers have focused their attention in their development.A relevant example are the works of Tadeu and Kausel 12 and of Tadeu and Ant ónio 13 , who proposed 2.5D Greens's functions for acoustic and elastic wave propagation in layered media, built as a superposition of the effects of plane waves with different inclinations; these functions have, in fact, been extensively used in subsequent works.Ant ónio et al. 14 developed a boundary element formulation incorporating Green's functions to describe 2.5 D wave propagation for the case of a waveguide with an elastic bottom and used them to study the scattering of waves by a buried or submerged object.A recent work by Pereira et al. 15 described a formulation based on the BEM which allows simulating the scattering of sound in an underwater configuration including a fluid seabed with multiple layers and a bottom discontinuity.
Recently, meshless methods have also been used for the study of underwater sound scattering in different types of environments.Different meshless methods are described in the literature, including Meshless-Local-Petrov-Galerkin MLPG methods see, e.g., Atluri 16 , RBF collocation methods see, e.g., Kansa 17,18 , or the method of fundamental solutions MFS e.g., Golberg and Chen 19 .Examples of the application of these strategies to underwater acoustics can be found in recent works by Godinho et al. 20 , using RBF base local interpolation methods formulated in the time domain, or by Costa et al. 21 , making use of the MFS together with the fundamental solutions for a flat waveguide and for a perfect wedge.
The scattering by a submerged object located within a fluid medium has also been investigated by researchers, and works describing the scattering features of submerged circular cylindrical elastic shell structures have also been published.The wave scattering by submerged elastic circular cylindrical shells, filled with air, struck by plane harmonic acoustic waves was analyzed by Veksler et al. 22 .In that work, the standard resonance scattering theory was used to study the modal resonances, focusing on the generation of bending waves.More recently, Godinho et al. 23 described an analytical solution for the scattering of such structures buried in a homogeneous fluid medium.Later, the same authors 24 used a BEM formulation to analyze the effect of a construction defect in the vibration of such structures.However, it is important to note that this BEM formulation degenerates whenever the thickness of the structure is very small, and therefore, alternative methods should be used.
In the present work, the authors address the case in which a regular circular shell structure is buried within a fluid seabed under a water-filled flat waveguide.The approach proposed here is based on a hybrid approach which incorporates the analytical solutions described in 23 for the submerged circular shell structures, together with the analytical solution known for a waveguide with a fluid bottom using the methodology proposed in 13 .The coupling of these solutions is performed in the fluid medium that describes the bottom by using the MFS and defining a virtual coupling boundary around the shell structure, along which the continuity of pressures and normal displacements is imposed.This formulation can easily incorporate multiple scattering objects, with different properties, although they are restricted to have a circular shape.More importantly, the method allows accounting for the full solid-fluid and fluid-fluid interaction that occurs at the physical interfaces of the system in an accurate manner, leading to precise results, since it is based on analytical solutions of each individual problem.Additionally, since it uses the analytical solution for a submerged circular shell, it allows modeling thin structures, overcoming the difficulties identified above for the BEM.
The paper is structured as follows: first, the governing equations of the problem are described in the frequency domain; then, the frequency domain multiregion MFS strategy for the coupling of the waveguide with the solid shells is formulated; there follows a description of the analytical solutions to be used for the submerged shell structures and for the waveguide with a fluid bottom; then, the proposed model is verified against BEM models; a procedure for obtaining time responses from the computed frequency-domain results is then described; finally, a numerical simulation is presented, illustrating the applicability of the model to a realistic configuration.

Governing PDEs
Within the scope of this work, the 2D scattering of waves by cylindrical shell structures embedded within a fluid medium is analyzed.Thus, the governing equations of the problem correspond to the vectorial and scalar wave equations, respectively, for the solid and for the fluid regions of the analysis domain.
Considering a homogeneous, linear isotropic elastic domain with mass density ρ s , shear wave velocity β s , and compressional wave velocity α s , the propagation of elastic waves can be described by vectorial wave equation where the vector u represents the displacement, ω is the circular frequency, and, for a twodimensional problem, ∇ ∂/∂x i ∂/∂y j; i and j are the unit vectors along the x and y directions.
If the propagation medium is a fluid, with mass density ρ f , the propagation is governed by the Helmholtz equation, which can be written as where p is the pressure and k f ω/α f is the wave number, with α f being the speed of sound in the fluid medium; for this scalar equation, ∇ 2 ∂ 2 /∂x 2 ∂ 2 /∂y 2 .Within this fluid medium, the displacements can be defined as a function of the first spatial derivative of p, and are given by 2.3

Formulation of the Hybrid Numerical-Analytical Model
Consider a fluid waveguide, with a fluid bottom simulating a sedimentary seabed.Within this seabed, consider the presence of an arbitrary number of circular, shell structures, made of elastic materials, and filled with a fluid material.This configuration is depicted in Figure 1 a .A hybrid analytical-numerical model based on the method of fundamental solutions MFS is proposed in this paper to calculate the pressure field within the waveguide generated by an acoustic source in the presence of such configurations.For this purpose, consider that in the presence of NR shell structures, the problem is divided in NR 1 subregions, one of them being the outer region, incorporating both the waveguide and the fluid bottom, and each of the NR subregions is defined around the shell structure, as represented in Figure 1 b .
If the fundamental solutions are known for each of the defined subregions, it becomes possible to establish a coupled model, which accounts for the full interaction between the involved fluids and the solids that compose the shell structures, by just establishing the continuity of pressures and displacements along the boundaries connecting the subregions.Using the MFS, the acoustic field in the outer subregion, containing the waveguide, can be defined by considering a number of virtual sources, NR j 1 NVS j placed within the remaining subregions, and combining their effects in a linear manner as a j,l G waveguide x, x vs j,l G waveguide x, x 0 , 3.1a while for a receiver placed within fluid of the jth inner subregion, we have where x represents a point of coordinates x, y , x 0 is the position of the real source illuminating the system, x vs j,l is the position of each of the NVS j virtual sources placed within subregion j, G waveguide x, x 0 is the fundamental solution for the waveguide with fluid bottom at a point x originated by a source positioned at x 0 ; G shell x, x 0 is the fundamental solution for each interior subregion, incorporating the full interaction between the shell structures and the outer and inner fluids; the coefficients a j,l are, "a-priori", unknown and must be determined by establishing a system of equations, enforcing the continuity of pressures and displacements along each of the NR boundaries separating the outer subregion from each internal subregion.Assuming that the boundary conditions are enforced at NVS j collocation points along the kth boundary as illustrated in Figure 1 b , the continuity equations along the mth collocation point x c,k m of that boundary can be written as where the coefficients b k,l are, "a-priori", unknown amplitudes of the fundamental solution for the region containing the shell structure.A N×N linear system of equations, with N 2× NR j 1 NVS j , can then be built.Once this system of equations is solved, one may obtain the pressure at any internal point by applying equations 3.1a and 3.1b .
An important point that should be highlighted concerning this formulation is that the coupling between subregions is enforced in fluid-fluid interfaces at some distance from the interfaces with the solid media that constitutes the shell structures.This strategy allows the coupling to be performed in a region with smooth variations of the pressure, which greatly improves the performance of the MFS.Additionally, since the interface between subregions is virtual, it can assume a smooth shape, such as that of a circle, which has been demonstrated in previous works that leads to very accurate results 25 .Finally, if the fundamental solutions are computed analytically within each subregion, a further step can be given towards obtaining high accuracy.In what follows, these fundamental solutions are described.

Analytical Solution for a Fluid Waveguide with a Fluid Bottom
Green's function for a flat fluid waveguide bounded bellow by a fluid halfspace simulating a seabed and above by a free surface can be obtained using the definition of displacement potentials, using the decomposition of the wavefield in terms of plane waves.These solutions are known for layered systems and can be derived following the methodology presented by Tadeu et al. 12, 13 .In this technique, the solutions can be expressed as the sum of the source terms equal to those in full space and of surface terms generated at the free surface and at the interface between the waveguide and the fluid halfspace.The calculation of the surface terms requires knowledge of the potentials' amplitudes.For the definition of these functions, consider the geometry depicted in Figure 2.
For an infinite fluid space the wavefield produced by a linear pressure load is here defined by a pressure potential as a superposition of plane waves by means of a discrete wavenumber representation, obtained by applying a Fourier transform along the x direction.The integrals of the expressions are transformed into a summation, by assuming an infinite number of virtual plane sources distributed along the x direction, at equal intervals, L x .To avoid the influence of neighboring fictitious sources in the response, complex frequencies of the form ω c ω − i × ξ are used here following the methodology described, for example, in 13 .This procedure allows to obtain the following pressure potential: where By adequate derivation of this potential, one obtains Green's function at point x for a infinite homogeneous fluid medium, when a pressure load is applied at x 0 , with coordinates x 0 , y 0 , as follows: The scattered wavefield in the waveguide can be defined in a similar way, by means of two displacement potentials, one representing the contribution of the top free surface φ 1 and the other related to the interface with the bottom halfspace φ 2 .These potentials are written as

Mathematical Problems in Engineering
To define the wavefield in the bottom halfspace, an additional potential must be defined in a similar manner, and is given by In these expressions B 1 n , B 2 n , and B 3 n are unknown coefficients to be determined after solving a system of equations, built so that the field, produced simultaneously by the source and the surface terms, should give the appropriate boundary conditions at the interfaces.The imposition of null pressure at the free surface, and of continuity of pressure and normal displacements at the fluid-fluid interface for each value of n, yields a system of three equations in the three unknowns see the appendix .Once the unknown coefficients have been calculated, the scattered pressure associated with the surface terms can be obtained.Green's functions for the fluid layer are then given by the sum of the source terms and the surface terms originated in both interfaces.
If a source acts in the top fluid waveguide , this leads to the following expressions for the pressure field in the system:

3.7
If the source is positioned in the seabed, a similar procedure can be used, including the source term in the pressure field of the bottom halfspace.From these equations, it becomes straightforward to apply 2.3 in order to obtain the displacements at any field point.

Analytical Solution for a Circular Cylindrical Shell Embedded in a Fluid Medium
Consider a circular shell solid structure, defined by the internal and external radii r A and r B , respectively, and submerged in a homogenous fluid medium, as illustrated in Figure 3.A harmonic dilatational source, placed in the exterior fluid medium, is assumed to illuminate the system, generating waves that hit the surface of the submerged structure.Part of the incident energy is then reflected back to the exterior fluid medium, with the rest being transmitted into the solid material, where they propagate as body and guided waves.These waves continue to propagate and may eventually hit the inner surface of the structure, where similar phenomena occur.The wavefield generated in the exterior fluid medium Fluid 2 depends both on the incident pressure waves and on those coming from the external surface of the shell.The latter propagate away from the cylindrical shell and can be defined using the following displacement potential when a cylindrical coordinate system is centered on the axis of the circular cylindrical shell, Inside the solid material of the shell, two distinct groups of waves exist, corresponding to inward travelling waves, generated at the external surface, and to outward travelling waves, generated at the internal surface of the shell.Each of these groups of waves can be represented using one dilatational and one shear potential where k αs ω/α s , k βs ω/β s , and α s and β s are, respectively, the dilatational and shear wave velocities permitted in the solid material J n • • • correspond to Bessel functions of the first kind and order n.
In the fluid that fills the shell structure Fluid 3 , only inward propagating waves are generated.For this case, the relevant dilatational potential is given by where k α3 ω/α f3 and α f3 is the pressure wave velocity in the inner fluid.The terms A j n j 1, 6 for each potential of 3.8 , 3.9 , and 3.10 are unknown coefficients to be determined by imposing the required boundary conditions.For this case, the necessary boundary conditions Mathematical Problems in Engineering are the continuity of normal displacements and stresses, and null tangential stresses on the two solid-fluid interfaces.
Consider, now, that the incident field, in terms of displacement potential, generated by the acoustic source located at x 0 , y 0 can be defined at a point x, y as

3.11
In order to establish the appropriate equation system, this incident field must also be expressed in terms of waves centered on the axis of the circular cylindrical shell structure.This can be achieved with the aid of Graf's addition theorem, leading to the expression in cylindrical coordinates where r 0 is the distance from the source to the axis of the circular cylindrical shell and ε n is 1 if n 0 and 2 in the remaining cases.
The solution of the equation system can then be used to compute the stresses in the solid medium as a summation of solutions obtained for pairs of values of n and k z .The final equation system can be found in published works, namely, 23, 24 .
After the solution of the corresponding equation system is computed, the unknown values A j n j 1, 6 can be used to determine the final wavefields.For the outer fluid, the pressure field at a point x, y can be written as

3.13
The corresponding displacement field can then be easily determined by applying 2.3 .

Calculation of Responses in the Time Domain
The pressure field in the spatial-temporal domain is obtained by modeling a Ricker wavelet, whose Fourier transform is in which Ω ωt o /2, A is the amplitude, t s is the time when the maximum occurs, and πt o is the characteristic dominant period of the wavelet.This wavelet form has been chosen, because it decays rapidly, both in time and frequency, reducing computational effort and allowing easier interpretation of the computed time series and synthetic waveforms.
The analysis uses complex frequencies, where ω c ω − iζ, with ζ 0.7Δω, which further reduces the influence of the neighboring fictitious sources and avoids the aliasing phenomena.In the time domain, this shift is later taken into account by applying an exponential window e ξt to the response 26 .

Model Verification
To verify the proposed coupling strategy, its results were compared with those obtained using alternative methodologies for a number of situations.Since no analytical solution is known for the complete problem to be solved here, this verification was performed against other numerical methods, namely, the boundary element method BEM .In what follows, two verification examples are described for specific cases, and a brief note on the stability of the procedure is presented.

Verification Example 1: A Circular Shell in an Halfspace Fluid Medium
In a first verification example, consider the case of an acoustic water halfspace, allowing a propagation velocity of 1500 m/s, and exhibiting a density of 1000 kg/m 3 , hosting a circular shell structure, made of an elastic material with a density of 1400 kg/m 3 , and allowing propagation velocities for the P and S waves of 2182.2 m/s and of 1336.6 m/s, respectively.This structure has an external radius of 1.5 m and an internal radius of 0.75 m, is filled with water, and is positioned with its centre at coordinates x 3.0 m and y −4.0 m.An acoustic source, located at x 0.0 m and y 5.0 m, illuminates this system, as illustrated in Figure 4 a .
The described configuration has been modeled using two different approaches.In the first, the proposed coupling strategy, making use of the fundamental solutions described above, is used.To allow the use of the fundamental solution for a waveguide over a fluid seabed, a virtual interface is considered at y 0.0 m, and the same properties are ascribed to the waveguide and to the fluid bottom.A virtual circular interface is also considered around the shell structure, with a radius of 1.6 m, in order to allow coupling the two fundamental solutions.80 collocation points are placed along this boundary, and two sets of 80 virtual sources are positioned as described in Section 3, at a distance of 0.5 m from the interface.The second model makes use of the boundary element method, as described in 24 , and a total of 240 elements is used to discretize the structure 160 for the outer boundary and 80 for the inner boundary .To account for the free surface of the halfspace, proper Green's functions are used for the outer medium.These functions are defined making use of the image source technique, where a single source is placed symmetrically to the real source with respect to the free surface, and with inverted polarity.
In both cases, the response is computed at a receiver positioned at x 6.0 m and y −3.0 m, for frequencies ranging from 20 Hz to 2500 Hz, with an increment of 20 Hz.Complex frequencies defined as, ω c ω − i × 0.7 × Δω are used in the calculation.
The results computed making use of both methods are depicted in Figure 4 b .As can be seen in this figure, the two sets of results match perfectly along the full set of frequencies analyzed.

Verification Example 2: Two Rigid Circular Inclusions Buried in a Fluid Seabed under a Fluid Waveguide
A second verification example has been analysed in order to assess the correctness of the results in the presence of more than two buried inclusions, positioned within a seabed with different properties from the waveguide.For this test, consider the case of an acoustic water waveguide medium properties are given in the previous section , with a depth of 20.0 m; bellow this waveguide, a fluid seabed is considered, allowing sound to travel at 2100 m/s, and exhibiting a density of 1800 kg/m 3 .Within the seabed, two circular rigid inclusions, with a radius of 0.5 m, are modeled.These inclusions are positioned at x 3.0 m and y −4.0 m and at x 6.0 m and y −4.0 m.An acoustic source, located at x 0.0 m and y 5.0 m, illuminates this system, as illustrated in Figure 5 a .

Mathematical Problems in Engineering
The described configuration has been modeled using two different approaches.Once again, the first corresponds to the proposed coupling strategy.Two virtual circular interfaces are considered around the circular inclusions, with a radius of 0.75 m, in order to allow coupling the fundamental solutions for the inner and outer domains.One should note that in this case, the fundamental solution for the case of a rigid inclusion can easily be obtained from Section 3.2, just considering the potential corresponding to the outer fluid and imposing the necessary null normal displacements at the outer interface.30 collocation points are placed along this virtual interface, and two sets of 30 virtual sources are positioned as described in Section 3, at a distance of 0.5 m from the interface.In a second model, the boundary element method is used, discretizing each of the two inclusions using 30 elements and the interface between the waveguide and the seabed using 950 elements.In order to limit the number of boundary elements used to discretize this interface, complex frequencies with an imaginary part are used ζ 0.7 2π/T .This considerably attenuates the contribution of the responses from the boundary elements placed at L 2α f T , reducing the length of the interface to be discretized see, e.g., 15 .In our calculations a value of T 0.05 s and L 210 m were used to define this discretization.The free surface of the halfspace is accounted using Green's functions defined by the image source technique.
In both cases, the response is computed at a receiver positioned at x 6.0 m and y −3.0 m, for frequencies ranging from 20 Hz to 1000 Hz, with an increment of 20 Hz.As  in the previous case, the results calculated making use of the two strategies are displayed in Figure 5 b .Again, the results match perfectly along the full set of frequencies analyzed here.

Behavior of the Coupled Model in the Presence of Two Circular Shell Structures Buried in a Fluid Seabed
An additional study was performed to better understand the behavior of the proposed model concerning the variability of its results with the number of collocation points and with the position of the virtual sources.For this purpose, consider the example illustrated in Figure 6, in which two buried elastic shell structures are embedded within a seabed with different properties from the waveguide.The properties of the acoustic water waveguide, 20.0 m deep, and of the fluid seabed are similar to those used in verification example described in the Section 5.2.The two circular structures have an external radius of 1.0 m and an internal radius of 0.5 m and are positioned at x 3.0 m and y −4.0 m and at x 6.0 m and y −4.0 m.The elastic properties of the shell structure are defined in Section 5.1.To couple the waveguide with the two structures, two virtual interfaces with a radius of 1.2 m are defined.The response has been calculated for a frequency of 2000 Hz, at a receiver placed at x 5.0 m and y −2.0 m, using different numbers of collocation points and positioning the virtual sources at different distances from the virtual interfaces between the shell regions and the waveguide region.
Table 1 presents the results at that receiver calculated for 30 collocation points when the distance between the virtual sources and the interface D assumes different values.In  that table, the relation D/R is used to define the distance as a function of the radius of the coupling interface R .As can be observed in that table, the response is stable as long as the virtual sources are not very close to the interface.In fact, for that case, a singularity of the fundamental solution occurs very close to the boundary, degrading the quality of the result.When D/R is 0.3 or larger, the variation of the response is very small and indicates a good behavior of the coupling strategy.Table 2 presents the results computed when D/R 0.4, and for varying numbers of collocation points.Here, the response can be seen to stabilize above 30 collocation points, corresponding to a relation between the wavelength and the distance between collocation points of just 5.This relation is relatively small when compared with those required for BEM discretizations around 7 to 8 .

Numerical Application
In order to illustrate the applicability of the proposed numerical approach, consider now a fluid waveguide, 20.0 m deep, with a sedimentary seabed, as displayed in Figure 7. Assume that the fluid inside the waveguide is water, with a density ρ f1 1000.0 kg/m 3 and allowing a dilatational wave velocity α f1 1500.0 m/s.The sedimentary seabed is first modeled with a density ρ f2 1800.0 kg/m 3 and permits the propagation of dilatational waves with a velocity α f2 2100.0 m/s.A second scenario is also considered, in which the seabed assumes a density ρ f2 1500.0 kg/m 3 and allows a sound velocity of α f2 1600.0 m/s, which corresponds to approaching the properties of the seabed to those of the fluid medium.

Mathematical Problems in Engineering
Within the seabed, consider the presence of two similar circular shell structures with external and internal radii r B 1.0 m and r A 0.95 m, respectively, made of an elastic material with density ρ s 7850.0 kg/m 3 , and allowing a dilatational wave velocity α s 6009.0 m/s and a shear wave velocity β s 3212.0 m/s; these structures are filled with a fluid material with the same properties as water.The described scenario is excited by a cylindrical pressure source placed near the bottom of the waveguide at x −10.0 m and y 2.0 m.
To simulate this problem, two virtual circular interfaces with a radius of 1.2 m are defined to account for the two subdomains containing the buried structures.Each of those interfaces is defined using 35 collocation points, and two sets of virtual sources are positioned at a distance of 0.4 times the radius of the virtual interface.
Calculations are then performed over a frequency range between 20.0 Hz and 2560.0Hz, assuming a frequency step of 20.0 Hz; for the purpose of calculating time responses, the defined increment allows a total analysis time of T 50.0 ms.Time domain signals are computed by means of an inverse Fourier transform, using the methodology described earlier.
The pressure field in the waveguide was computed over a grid of receivers, equally spaced at Δx 0.5 m, Δy 0.5 m, placed between x 0.0 m and x 25.0 m and y 0.0 m and y 20.0 m.A sequence of snapshots displaying the pressure wave field over the grid of receivers, at different instants, is presented to illustrate the results.Figure 8 displays snapshots of the pressure response, for different time instants, over the grid of receivers placed in the waveguide, generated by a source emitting a Ricker pulse with a characteristic frequency f k 400 Hz.A grayscale is used to represent the amplitudes of the waves arriving at the receivers, with lighter colors corresponding to higher values and darker colors representing lower values.These responses were computed assuming the waveguide with a sedimentary seabed which allows the propagation of dilatational waves with a velocity α f2 2100.0 m/s without shell structures buried.
At time t 0.0 ms, the load creates a cylindrical pressure wave that propagates away from it.In the snapshot of Figure 8 a , corresponding to t 14.6 ms, this incident pulse is visible identified as P1 , followed by a first reflection from the bottom of the waveguide identified as P2 .At receivers placed near the ground, a third reflection may also be identified, which is related to the head wave generated in the surface of the seabed identified as P3 .This wave is originated at the interface between the two media and travels along this interface with the velocity of the faster medium, which is the seabed, with α f2 2100.0 m/s; therefore, it appears in the plot at receivers placed farther from the source.As time increases, it is possible to identify the reflections generated at the free surface identified as P4 , with inverted polarity see Figure 8 b .For subsequent instants, a sequence of pulses originated by multiple reflections in the surface and bottom of the waveguide can be identified see Figures 8 c and 8 d .These reflections tend to lose energy as time increases, with part of the energy being transmitted to the seabed, and a stationary field is generated inside the waveguide by these waves, which travel up and down between the surfaces of the channel and tend to become flat as time increases.
When the two shell structures are buried in the seabed a different wave pattern inside the waveguide may be created.In order to assess the presence of these structures under different conditions, snapshots of the sound propagation within the waveguide were captured for two different sets of properties of the seabed: α f2 2100.0 m/s see first column of Figure 9 and α f2 1600.0 m/s see second column of Figure 9 , respectively.When the seabed allows a velocity α f2 2100.0 m/s see first column of Figure 9 , a set of additional pulses appear in the response labeled as Pshell , which refer to reflections originated by the presence of the shell structures.These reflections can further be identified in the snapshots corresponding to subsequent instants see Figures 9 b1 -9 d1 although displaying smaller amplitudes, due to the contrast between media, which tends to hinder energy exchanges.As expected, when the seabed assumes a dilatational wave velocity which approaches that of the fluid in the waveguide see plots provided in the second column of Figure 9 , the amplitudes of the scattered pulses provided by the shell structures are increased, providing a clear perception of their presence.For later instants, the responses display multiple pulses, related not only to reflections of waves generated in the waveguide, but also to several reflections originated at the shell structures, at the top and at the bottom of the waveguide.It is also interesting to note that for both cases, the reflection pattern originated at the buried structures is quite complex, revealing multiple reverberation effects that occur not only between the structures and the sea bottom, but also within the structures themselves and within the fluid that fills their interior.
With the aim of understanding how the properties of the fluid inside the shell structures may influence the responses, the time domain response, originated by a source emitting a Ricker pulse with a characteristic frequency f k 800 Hz, has been calculated along a line velocity of α f2 2100.0 m/s.Responses computed for the waveguide without the buried structures are displayed in Figure 10 a as a reference.In that first plot, the initial reflections occurring between the free surface and the fluid-fluid interface can clearly be identified, with a 180 • phase change occurring when the incident pulse is reflected at the free surface P2 and P4 .The so-called "head wave" can also be identified in this plot P3 , travelling at 2100 m/s along the interface.When water-filled shell structures are introduced Figure 10 b a clear identification of a sequence of scattered pulses originated by these structures is also possible.Interestingly, this sequence allows identifying a reverberation effect in which multiple reflections are occurring within the structure, both on the solid shell and on the filling fluid.This effect is even clearer in the second series of pulses, originated when the pulses coming from the free surface hit the buried structure.In fact, as the incident pulse hits the interface at an almost tangent angle, it is mainly reflected back to the waveguide, and very little energy penetrates the seabed; this no longer occurs for the second set of pulses, which reach this surface at larger angles, and thus allow more energy to be transmitted to the bottom.Some differences are visible when the shell structure is filled with air, for which case the scattered pulses are captured at these receivers with smaller amplitudes; in fact, the high contrast between the solid and the interior fluid air , makes the energy exchanges among materials more difficult to happen.As a consequence, energy tends to dissipate faster, as time increases, and thus reduces significantly the reverberation effect.In the reference result of Figure 10 d , considering the external boundary of the inclusions to be rigid, no energy propagates to the interior of the shell, and, as a result, the response reveals fewer pulses coming from the buried structures.These results clearly show the importance of accurately modeling both the solid of the shell and the fluid in its interior and show that the commonly used approximation of assuming a rigid behavior of the structures neglects important parts of the propagation phenomena.
To emphasize these findings, two further plots are presented in Figure 11, considering the water-filled and air-filled structures buried in a seabed which allows a propagation velocity of 1600 m/s.In Figures 11 a and 11 b , plots are presented illustrating the time responses for water-filled shells and air-filled shells.For these cases, the difference between the scattering patterns is significantly increased, due to the smaller contrast between the waveguide and the seabed.The reverberation effect is now very clear in the case of waterfilled shells, with rings of pulses being generated due to multiple reflections within the structure and filling fluid.This effect is much less pronounced when the structure is filled with air.One additional feature of this response corresponds to the presence of two distinct sets of pulses, originated at each of the two inclusions; the first arrival of each set occurs approximately at the receivers placed immediately above the buried structure and allows an identification of the presence of the two separate shells.

Conclusions
In this paper, the coupling between different analytical solutions using the MFS is proposed to address the problem of scattering of acoustic waves in a waveguide in the presence of buried structures.The scattering structure is assumed to be buried in the fluid seabed bellow a water waveguide and is a circular elastic shell filled with a fluid that may have different properties from the host medium.The proposed strategy was formulated and implemented and was shown to provide good results when compared with alternative numerical modeling techniques.Since it performs the coupling between closed-form solutions, the method provides accurate results, while allowing a compact and simple model description.One major advantage of the proposed model is that it allows the simulation of very thin solid structures, without the problems usually associated with thin bodies when using alternative methods.A number of applications were presented, revealing that taking into account both the elastic properties of the buried shell and the properties of a fluid which fills that structure can be important, leading to marked differences when results are compared with usual simplifications, such as considering the buried structure to be rigid.

Appendix
This appendix presents the system of equations required to obtain Green's function for a flat fluid waveguide bounded bellow by a fluid halfspace and above by a free surface as defined in Figure 2.
This system can be defined as

Figure 1 :
Figure 1: a Schematic representation of the problem to be solved; b Detail of the interface around one of the shell structures, indicating the distribution of virtual sources and the coupling boundary interface.

Figure 2 :
Figure 2: System with a fluid waveguide over a fluid seabed.

− x 0 2 y
− y 0 2 and H n • • • represent Hankel functions of the second kind and order n.

Figure 3 :
Figure 3: Circular cylindrical shell structure submerged in a fluid medium.

H = 20 = 1000 kg/m 3 αFigure 4 :
Figure 4: a Schematic representation of the first verification example.b Responses provided by the BEM and by the proposed coupled numerical-analytical model.

H = 20 Figure 5 :
Figure 5: a Schematic representation of the second verification example.b Results provided by the BEM and by the proposed coupled numerical-analytical model.

H = 20 Figure 6 :
Figure 6: Schematic representation of the system with two buried shell structures.

H = 20 Figure 7 :
Figure 7: Geometry defined for the numerical applications.

Figure 8 :
Figure 8: Snapshots displaying the pressure wave field over the grid of receivers at different instants, in the waveguide assuming the seabed with α f2 2100.0 m/s, without the shell structures: a t 14.6 ms; b t 24.4 ms; c t 34.2 ms; d t 39.1 ms.

Figure 9 :Figure 10 :
Figure 9: Snapshots displaying the pressure wave field over the grid of receivers placed in the waveguide with a seabed, where shell structures are buried: a t 14.6 ms; b t 24.4 ms; c t 34.2 ms; d t 39.1 ms.

Figure 11 :
Figure 11: Time domain responses, captured along a line of receivers placed at y 0.5 m, assuming the waveguide with a seabed, allowing a dilatational wave velocity of α f2 1600.0 m/s and with a the shell structures filled with water and b the shell structures filled with air.

Table 1 :
Response at the field point x 5 m; y −2 m when 30 collocation points are used and for different positions of the virtual sources.

Table 2 :
Response at the field point x 5 m; y −2 m when different numbers of collocation points are used and the distance from the virtual sources to the interface is 0.4 times the radius of the fictitious interface.