Figure 1 depicts a system involving the interaction between acoustic and structural components. An elastic thin-shell occupies the region Ω s , which is encompassed by an infinite fluid domain Ω f . The interface where the fluid and the structure meet is denoted by Γ s f , defined as the intersection of Ω s and Ω f . The space within the shell is occupied by air, forming the domain Ω a . Given the significant difference in density between the air and the structural material, it is reasonable to disregard the influence of air pressure on the structural vibrations.

The mechanical response of the thin-shell is defined by the Kirchhoff-Love shell theory, as detailed in references [36, 37]. Meanwhile, the acoustic fields within the fluid domain are regulated by the Helmholtz equation, which is discussed in [38-41]. The thin-shell structure is actuated by a time-varying force with an angular frequency of ω . Consequently, the equations governing the acoustic-structural interaction system are presented as Equations 1–4. ∇ ⋅ σ x + ω 2 ρ s u x = 0 x ∈ Ω s (1) ∇ 2 p x + k 2 p x = 0 x ∈ Ω f (2) σ x ⋅ n s x − p x n f x = 0 x ∈ Γ s f (3) − i ω u x ⋅ n f x = v f n x x ∈ Γ s f (4)

Equations 1, 2 serve as the foundational equations that shape the configuration of acoustic fields. Equations 3, 4 delineate the conditions for continuous displacement and balanced forces at the interfaces where the fluid and structure interact. The symbol ∇ denotes the Laplacian operator, σ signifies the stress tensor, u represents the displacement vector, ρ s is the density of the shell, p ( x ) denotes the acoustic pressure, n s is the outward normal vector to the shell surfaces, n f is the inward normal vector to the structural surface, i is the imaginary unit with i = − 1 , v f n is the normal velocity of the fluid, and k , defined as ω / c , represents the wave number, with c being the speed of sound propagation within the fluid domain.

The choice of a discretization approach or the selection of basis functions is crucial in numerical simulations. Typically, in traditional FEM and BEM, polynomial functions are the preferred choice for basis functions. In this study, we integrate FEM and BEM within the framework of isogeometric analysis (IGA), as referenced in [42-46]. This integration leverages the same basis functions used for geometric modeling to discretize the physical fields. This approach not only preserves the geometric precision but also eliminates the need for mesh generation. Most significantly, it facilitates the construction of high-order continuous fields, which is essential for the application of the Kirchhoff-Love shell theory, as further elaborated in [47].

Non-Uniform Rational B-splines (NURBS), as discussed in [48- 52], are a prevalent choice for geometric construction in Computer-Aided Design (CAD), with polygonal meshes subsequently derived for numerical analysis. In this research, we embrace the principles of isogeometric analysis, utilizing identical basis functions for both geometric modeling and numerical simulation. We opt for Catmull-Clark subdivision surfaces over NURBS for CAD model construction, as they ensure watertight geometries with support for arbitrary topologies, as described in [53].

Constructing geometry with Catmull-Clark subdivision surfaces begins with the establishment of a control mesh, which is made up of quadrilateral elements. The vertices of this mesh are known as control points. Following the initial setup, the control mesh undergoes subdivision, during which new control points are introduced and the positions of the existing ones are adjusted, as shown in Figure 2. This process of subdivision can be repeated iteratively, adhering to a set of subdivision rules that progress from level k to level k + 1 , as illustrated in Figure 3. However, since continuous subdivision is not feasible, practical implementations of limit subdivision surfaces are assessed through the use of spline basis functions for parameterization.

Subdivision surfaces excel in their capacity to manage singular points while ensuring controlled curvature at these locations. Figure 4A illustrates a regularly colored subdivision surface element, where each vertex within the element has a valence of 4, signifying the lack of extraordinary vertices within the element. A patch is delineated as a collection of all elements that share vertices with a particular target element. A regular patch is composed of nine elements featuring sixteen vertices. The element highlighted in Figure 4B contains an extraordinary vertex, which adds complexity to surface evaluation as it disrupts the tensor-product property. It should be noted that the initial grid might include irregular elements with multiple extraordinary vertices, but after subdivision, each resulting sub-element will contain at most one exceptional vertex. The evaluation of a surface point with parametric coordinates ( ξ , η ) is conducted through a linear combination of control points and 2 κ + 8 basis functions, characterized by a valence of κ , as shown in Equation 5. x ξ , η = ∑ ℓ = 0 2 κ + 7 N ^ ℓ ξ , η P ℓ (5) in which x represents the Cartesian coordinates of the point, N ̂ ℓ refers to the basis functions, P ℓ signifies the control points, and κ indicates the valence of the vertex.

As previously stated, subdivision surfaces are employed for both crafting geometric models and discretizing the physical fields. The collective elements establish the boundary Γ , which is delineated by Equation 6. Γ = ⋃ e = 1 N e Γ e (6) in which N e denotes the overall count of elements, with element Γ e signifying the e -th element. Consequently, the displacement can be articulated as follows by leveraging the basis functions associated with the Catmull-Clark subdivision, as shown in Equation 7. u ̃ e = ∑ ℓ = 0 2 κ + 7 N ^ ℓ ξ , η u ̃ ℓ e (7) in which u ̃ ℓ e denotes the nodal parameter associated with the ℓ -th control point, employed for discretizing the displacements. Combining all the individual element matrices results in the subsequent Equation 8. A u ̃ = f (8) in which u ̃ represents the vector of displacements at all vertices, and f signifies the load vector. Furthermore, the matrix A can be formulated as Equation 9. A = K − ω 2 M , K = ⋃ e = 1 N e K e , M = ⋃ e = 1 N e M e (9) in which K signifies the global stiffness matrix of the structure, while M represents the global mass matrix. The element stiffness matrix is denoted by K e and the element mass matrix by M e . It should be noted that the load vector f is the sum of two components, f s and f p , where f s corresponds to the external mechanical load, and f p corresponds to the load due to acoustic pressure.

Within the three-dimensional domain Ω ⊂ R 3 bounded by Γ , the propagation of acoustic waves is governed by the Helmholtz equation. This equation can be transformed into boundary integral equations. To address the issue of irregular frequencies or spurious modes in exterior acoustic scenarios, we employ the Burton-Miller formulation, which is a linear combination of the conventional boundary integral equation (CBIE) and its derivatives [54-57], as shown in Equation 10. C x p x + α C x q x = ∫ Γ G x , y q y d Γ + α ∫ Γ ∂ G x , y ∂ n f x q y d Γ − ∫ Γ ∂ G x , y ∂ n f y p y d Γ − α ∫ Γ ∂ 2 G x , y ∂ n f x ∂ n f y p y d Γ + p ̃ inc x (10) in which x and y correspond to the source and field points, respectively. The term C ( x ) represents the discontinuity across the boundary, which takes a value of 1 2 for smooth surfaces. The coupling parameter α is given by i k when k > 1 , and α = i in other cases. The quantity q ( y ) = ∂ p ( y ) ∂ n f ( y ) denotes the flux of sound pressure. The term p ̃ inc ( x ) = p inc ( x ) + ∂ p inc ( x ) ∂ n f ( x ) encompasses both the incident wave’s sound pressure and its normal derivative at point x , with p inc ( x ) being the incident wave’s sound pressure at that location. The functions G ( x , y ) and its normal derivative for three-dimensional full-space acoustics are specified as Equations 11, 12. G x , y = e i k r 4 π r (11) ∂ G x , y ∂ n f y = − e i k r 4 π r 2 1 − i k r ∂ r ∂ n f y (12) in which r = | x − y | . The fundamental solution for half-space problems is represented by G * ( x , y ) , as shown in Equation 13. G * x , y = G x , y + β G x * , y (13) in which x * denotes the reflection of point x across the infinite plane Γ h , and β signifies the reflection coefficient. The coefficient β equals one for a rigid boundary and −1 for a soft boundary. By replacing G ( x , y ) with G * ( x , y ) , Equation 10 is converted into Equation 14. C x p x + α C x q x = ∫ Γ G * x , y q y d Γ + α ∫ Γ ∂ G * x , y ∂ n f x q y d Γ − ∫ Γ ∂ G * x , y ∂ n f y p y d Γ − α ∫ Γ ∂ 2 G * x , y ∂ n f x ∂ n f y p y d Γ + p ̃ inc * x (14) in which p ̃ inc * ( x ) = p inc * ( x ) + ∂ p inc * ( x ) ∂ n f ( x ) , and the partial derivative of G*(x,y) is shown in Equations 15–17. ∂ G * x , y ∂ n f y = − e i k r 4 π r 2 1 − i k r ∂ r ∂ n f y − β e i k r * 4 π r * 2 1 − i k r * ∂ r * ∂ n f y (15) ∂ G * x , y ∂ n f x = − e i k r 4 π r 2 1 − i k r ∂ r ∂ n f x − β e i k r * 4 π r * 2 1 − i k r * ∂ r * ∂ n f x (16) ∂ 2 G * x , y ∂ n f x ∂ n f y = e i k r 4 π r 3 3 − 3 i k r − k 2 r 2 ∂ r ∂ n f y ∂ r ∂ n f x + 1 − i k r n i x n i y + β e i k r * 4 π r * 3 3 − 3 i k r * − k 2 r * 2 ∂ r ∂ n f y ∂ r ∂ n f x + 1 − i k r * n i x n i y (17) In this approach, the physical field is discretized using the basis functions of the subdivision surface that are employed in the geometric construction, as displayed in Equation 18. p e = ∑ ℓ = 0 2 κ + 7 N ^ ℓ ξ , η p ̃ ℓ e q e = ∑ ℓ = 0 2 κ + 7 N ^ ℓ ξ , η q ̃ ℓ e (18) in which p e represents the acoustic pressure at the point ( ξ , η ) within the element Γ e ; q e signifies the normal flux associated with p e . Additionally, p ̃ ℓ e refers to the ℓ -th nodal variable associated with the sound pressure, while q ̃ ℓ e represents the normal flux corresponding to p ̃ ℓ e . By inserting Equation 18 into Equation 14 and applying these equations at a discrete set of collocation points x j , Equation 19 can be derived. C x j p x j + α C x j q x j + ∑ e = 1 N e ∑ ℓ = 0 2 κ + 7 ∫ Γ e N ^ ℓ ξ , η ∂ G * x j , y ξ , η ∂ n f y ξ , η d Γ p ̃ ℓ e + α ∑ e = 1 N e ∑ ℓ = 0 2 κ + 7 ∫ Γ e N ^ ℓ ξ , η ∂ 2 G * x j , y ξ , η ∂ n f x j ∂ n f y ξ , η d Γ p ̃ ℓ e = ∑ e = 1 N e ∑ ℓ = 0 2 κ + 7 ∫ Γ e N ^ ℓ ξ , η G * x j , y ξ , η d Γ q ̃ ℓ e + α ∑ e = 1 N e ∑ ℓ = 0 2 κ + 7 ∫ Γ e N ^ ℓ ξ , η ∂ G * x j , y ξ , η ∂ n f x j d Γ q ̃ ℓ e + p ̃ inc * x j (19) in which the subscript j indicates the index of the collocation point. These collocation points are determined by projecting the control points onto the surface. Equation 19 can be expressed in matrix-vector form as Equation 20. H p ̃ = G q ̃ + p i (20) in whichn H and G are the matrices of coefficients, while p ̃ and q ̃ are the column vectors that compile the nodal parameters related to the acoustic pressure and its flux, respectively. p i signifies the nodal pressure vector attributed to the incident wave. It is crucial to address the singular integrals present in Equation 10. Numerous methods have been proposed to handle these singular integrals within BEM, as cited in [58-61]. This research adopts the singularity subtraction technique introduced by Guiggiani, recognized for its precision and efficiency, as detailed in [62].

Equations 8, 20, originating from the structure (using FEM) and the acoustic (using BEM) respectively, are tightly interlinked and cannot be resolved in isolation. They are interconnected through the boundary conditions outlined in Equations 3, 4. The sound pressure within the fluid domain can be perceived as a force exerted on the shell surface. Consequently, the nodal force vector f p due to the acoustic load can be formulated as Equation 21. f p = n ̃ f Θ p ̃ (21) in which Θ is defined as the integral over the boundary Γ of the product of the transpose of the subdivision basis functions’ global vector N ^ ℓ and N ^ ℓ itself, multiplied by d Γ . The vector N ^ ℓ denotes the global vector of subdivision basis functions, and n ̃ f represents the matrix of normal vectors, which is given by Equation 22. n ̃ f = n 1 f ⋅ e 1 0  … n 1 f ⋅ e 2 0  … n 1 f ⋅ e 3 0  … 0  n 2 f ⋅ e 1 … 0  n 2 f ⋅ e 2 … 0  n 2 f ⋅ e 3 … … … … (22) in which e 1 , e 2 , e 3 denote the three unit normal vectors. By applying Equation 21, we can derive a new expression for the global nodal force vector, which is presented in Equation 23. f = C s f p ̃ + f s (23) in which C s f = n ̃ f Θ .

Next, we will investigate the relationship between the velocities on the shell’s mid-surface and the sound pressure. The vector v f n represents the nodal parameters related to the fluid’s normal velocity component, while the vector v s n signifies the nodal variables associated with the normal velocity component of the structure. Assuming no energy dissipation occurs at the interface between the fluid and structural domains, we get Equation 24. v f n − v s n = 0 (24) in which v s n = i ω C f s u ̃ , and C f s = n ̃ f T . Subsequently, the vector q ̃ , encompassing the nodal parameters related to the acoustic flux, can be formulated as Equation 25. q ̃ = ∂ p ̃ ∂ n = − i ω ρ f v f n = ω 2 ρ f C f s u ̃ (25) in which ρ f represents the fluid’s density.

By inserting Equation 25 into Equation 20, we obtain the following coupled system of equations for the acoustics, as shown in Equation 26. H p ̃ = G ω 2 ρ f C f s u ̃ + p i (26) Incorporating Equation 23 into Equation 8 yields the following coupled system of equations for structural dynamics shown in Equation 27. A u ̃ = C s f p ̃ + f s (27) By integrating Equation 27 into Equation 26, the indeterminate displacement u ̃ is eliminated, leading to a boundary element formulation that is linked with the finite element equation, as shown in Equation 28. H - G Y p ̃ = G q ̃ s + p i (28) in which Y denotes the global admittance matrix, and q ̃ s accounts for the influence of sound speeds originating from the structural domain. These are given by Equation 29. Y = ω 2 ρ f C f s A − 1 C s f q ̃ s = ω 2 ρ f C f s A − 1 f s (29)

This section examines the acoustic field scattered by a spherical shell structure when subjected to incident plane waves, taking into account the reflections from the seabed, as depicted in Figure 6. In the diagram, the center of the sphere is at the origin; Λ represents the seabed surface; r ̃ is the distance from the origin to the seabed surface; f p is the incident plane wave; Ω f signifies the infinite fluid domain. The seabed reflection coefficient is denoted by β , where β = − 1 corresponds to a soft seabed and β = 1 corresponds to a rigid seabed.The incident plane waves, which have an amplitude of 1, travel along the positive x -axis. Using Catmull-Clark subdivision surfaces, three elastic spherical shells with 3,458, 4,800, and 9,408 mesh points were constructed. Figure 7 presents the variation curves of sound pressure with frequency for the spherical models at 3,458, 4,800, and 9,408 mesh points, respectively. The computation times for a single frequency were 33 s, 56 s, and 85 s for each model. It can be observed from the figure that the differences in the sound pressure curves across the three mesh densities are minimal; however, the computation times differ significantly. Therefore, this study adopts the spherical model with fewer mesh points to balance computational efficiency and accuracy.When analyzing the acoustic scattering of the sphere model, we found that the computation time for the traditional Catmull-Clark subdivision surfaces coupled with the FEM-BEM method was 428 s, whereas the computation time for the Catmull-Clark subdivision surfaces accelerated by FMM and coupled with the FEM-BEM method was only 33 s, demonstrating that the computational efficiency of FMM is significantly higher than that of traditional algorithms.

The sound pressure values at the computed points on the circle for incident wave frequencies of 100 Hz and 200 Hz are depicted in Figure 8, where the circle represents the projection of the sphere model, with a radius of 2 m, onto the y o z plane. As observed in the figure, at a frequency of 200 Hz, the sound pressure curve exhibits more pronounced fluctuations, indicating that the larger the external sound wave impacting the sphere model, the greater the magnitude of sound pressure variation. When β = 0 , the sound pressure curves from numerical and analytical solutions align, confirming the validity and reliability of the algorithm presented in this study. For β = 0,1 , − 1 , the numerical solutions vary significantly, demonstrating the importance of considering the seabed reflection effects and highlighting the significance of the seabed’s nature.

Figure 9 illustrates the displacement and sound pressure distributions on the sphere model when it is positioned 5 m above the seabed. The first three columns represent the displacement components in the x , y , and z axes, respectively, while the last column indicates the sound pressure. As observed in the figure, the sound pressure pattern exhibits symmetry along the x -axis, and the displacement components also mirror this symmetry. The displacement in the x -direction is the most significant, whereas the y -direction shows the least. The sound pressure on the sphere’s surface diminishes as we move in the negative x -axis direction, whereas the x -direction displacement component increases progressively.

Figure 10 displays the sound pressure of the sphere model across various frequencies with β = 0 . Figures 10A–C depict the real part, imaginary part, and magnitude of the sound pressure, respectively. From the figure, it can be seen that when the frequency is within the range of [0.140]Hz, the analytical solution and numerical solution of the sound pressure on the surface of the sphere model match very well, and the sound pressure changes are very small. When the frequency is within the range of [160,200] Hz, the numerical solution still matches well with the analytical solution, but the sound pressure begins to change dramatically. This phenomenon further validates the reliability of the algorithm presented in this paper and shows that sound scattering is frequency-dependent.

Figure 11 shows the sound pressure of the sphere model at different frequencies when β = 0 , β = − 1 , and β = 1 . From the figure, it can be seen that the sound pressure curves for β = 0 , β = − 1 , and β = 1 are distinct, indicating that the impact of seabed reflection on the structural acoustic coupling system response in shallow water environments cannot be ignored. It is noteworthy that the sound pressure curves for β = 1 and β = − 1 are roughly distributed on either side of the sound pressure curve for β = 0 , and the two curves are essentially symmetric about the sound pressure curve for β = 0 . This phenomenon indicates that the effects of rigid and flexible seabeds on sound waves are opposite.

Figure 12 illustrates the sound pressure of the sphere model at various distances from the seabed, where r ̃ denotes the distance from the model’s center to the seabed. It can be observed from the figure that when the frequency is less than 60Hz, the sound pressure curves for different r ̃ are very close to each other; when the frequency exceeds 60Hz, the differences between the sound pressure curves for different r ̃ gradually increase. This phenomenon confirms that the effect of ground reflection is frequency-dependent.

This section analyzes the acoustic field scattering by a fish model under the action of incident plane waves, taking into account the reflection effects of the seabed, as shown in Figure 13. This subsection calculates the sound pressure and displacement under incident plane waves at different frequencies. The incident plane waves propagate along the positive x-axis with a unit amplitude and are scattered by the underwater model.

Using Catmull-Clark subdivision surfaces, two fish models with 54,216 and 222 mesh points were constructed. Figure 14 illustrates the sound pressure curves as a function of frequency for the fish models with 54,216 and 222 mesh points, respectively. The computation times for a single frequency were 483 s and 1,016 s for each model. As shown in the figure, the differences in the sound pressure curves between the two mesh densities are minor, with only slight fluctuations that are acceptable considering the significant disparity in computation times. Therefore, this study selects the fish model with 54,216 mesh points to optimize computational efficiency while maintaining acceptable accuracy.

When analyzing the acoustic scattering of the fish model, we found that the computation time for the traditional Catmull-Clark subdivision surfaces coupled with the FEM-BEM method was 5,143 s, whereas the computation time for the Catmull-Clark subdivision surfaces accelerated by FMM and coupled with the FEM-BEM method was only 483 s.

Figure 15 shows the distribution of displacement and sound pressure for the fish model at a distance of 0.5 m from the seabed. As shown in the figure, the sound pressure and displacement distributions are symmetric along the x-axis. The magnitude order of the various displacement components of the fish model at different frequencies is not the same. For instance, at a frequency of 1,000 Hz, the displacement component in the x-direction is the largest, and the displacement component in the z-direction is the smallest; at a frequency of 3,000 Hz, the displacement component in the z-direction is the largest, and the displacement component in the y-direction is the smallest.

Figure 16 displays the sound pressure of the fish model at different frequencies for β = 0 , β = − 1 , and β = 1 , with a frequency range of [03,000]Hz and a frequency step of 10 Hz. As shown in the figure, the sound pressure on the surface of the fish model varies with β . When the frequency is within the range of [0.600]Hz, the sound pressure curves for β = 0 , β = − 1 , and β = 1 match very closely. However, when the frequency exceeds 600Hz, differences begin to emerge among them. This phenomenon once again confirms that the seabed reflection effect is frequency-dependent, not only for spherical models but also for complex models, and it also validates the necessity of considering the seabed reflection effect.

Figure 17 illustrates the sound pressure of the fish model at various distances from the seabed. It can be observed from the figure that as the distance from the seabed increases, the sound pressure curve gradually aligns with the sound pressure curve for β = 0 . This phenomenon indicates that in shallow sea areas, the seabed reflection effect must be taken into account, while in deep sea areas, the seabed reflection effect can be selectively considered based on needs. This phenomenon is related to the nature of the seabed.