Since its introduction in the 1970s, subdivision surfaces have been extensively utilized in computer animation and graphics [54]. Additionally, they are accessible in the majority of industrial CAD solid modeling programs (e.g., CATIA, Creo). Subdivision surfaces are often seen in computer graphics literature as a method for repeatedly smoothing and refining a control mesh in order to generate smooth limit surfaces. They may also be seen from the perspective of finite and boundary element analysis as the generalization of splines to arbitrarily connected meshes.

The concept of subdivision schemes is to create a smooth surface out of a rough polygon mesh. Subdivision refinement systems, which may be categorized as interpolating or approximation schemes, create a smooth surface by a limited process of repetitive refinement beginning with an initial control mesh. All control meshes produced during subdivision refinement describe the exact same spline surface since the subdivision surfaces inherit the refinability attribute from the splines.

In this research, the structural-acoustic coupling analysis is conducted using a Loop subdivision technique [42].

The BEM and FEM are used to simulate the fluid and structural subdomains, respectively. The discretized boundary integral formulation of the fluid solution to the Helmholtz problem is given by [55], as shown in Equation 1: H p = G q + p i , (1) where

 H and G are the frequency-dependent BEM influencing matrices,

 p is the vector containing the nodal values for pressure,

 q is the vector of normal derivatives of p,

 p _i is the incident wave’s nodal pressure.

The continuous linear and quadratic element technique is often used, and alternatives to discontinuous elements with good accuracy have already been examined [14]. Discontinuous boundary elements outperform continuous boundary elements [14]. Interpolation nodes are positioned within discontinuous boundary elements, and the expressions of the interpolation functions depend on the position of the node within the element. Thus, by varying the placements of the interpolation nodes, a numerical solution with varied calculation accuracies may be achieved. Both discontinuous and continuous components are employed in the numerical computation in this study.

The steady-state response of the structure may be derived from the frequency-response analysis when a harmonic load with a transitory function e^−iωt is applied to it. When the acoustic-structure interaction is addressed, Equation 2 derives the linear system of equations to compute the nodal displacements u [51]. K − ω 2 M + i ω C u = f s + C sf p , (2) C sf = ∫ Γ N s T n N f d Γ , (3) where

 K is the structure’s stiffness matrix,

 M is the structure’s mass matrix,

 C is the structure’s damping matrix,

 ω is the angular frequency,

 f _s is the vector of the nodal structural forces,

 C _sf p is the acoustic load.

 Γ is the coupling boundary surface,

 N _s is the global interpolation function for the structural domain,

 N _f is the global interpolation function for the fluid domain,

 n is the normal vector of the surface.

The continuity condition across the interaction surface (Fritze et al.2005) is introduced in Equation 4, as follows: q = i ω ρ f v f , (4) v f = − i ω S − 1 C fs u , (5) S = ∫ Γ N f T N f d Γ , (6) C fs = C sf T , (7) where

ρ _f is the density of the fluid,

v _f is the normal velocity vector of the fluid.

Combing Equations 1, 2 and 4, Equation 8 gives the coupled system of equations of an elastic structure submerged in a heavy fluid derived. K − ω 2 M + i ω C − C sf − ω 2 ρ f G S − 1 C fs H u p = f s p i . (8)

On Equation 8, the application of an iterative solution (e.g., generalized minimum residual approach (GMRES)) results in unsatisfactory convergence. Substituting the finite element formulation into the boundary element equation to obtain a simplified system equation is an appropriate strategy [51]; [44, 56, 57], as shown in Equation 9: H p − G W C sf p = p i + G W f s , (9) W = ω 2 ρ f S − 1 C fs A − 1 , (10) A = K − ω 2 M + i ω C . (11)

It takes much time to solve A ⁻¹ directly in Equation 9. In fact, A sparse direct solver could be used to readily solve this symmetric, frequency-dependent system of equations Ax = f _s , thus obtaining the term A ⁻¹ f _s in Equation 9 quickly. The iterative solver GMRES [24] is introduced in this work to speed up the computing of solutions to the equations for the coupled boundary element system. There is no need to solve A ⁻¹ C _sf directly on the left-hand side of Equation 9. Considering the current iterative solution p _k for vector p in Equation 9, when a sparse direct solver is employed to solve the symmetric and frequency-dependent system of the linear equation Ax = C _sf p _k , the solution of A ⁻¹ C _sf p _k could be achieved effectively. Based on this, we could derive the solution of vector u by solving Eq. 9 and inserting the solution of vector p into Equation 2.

When dealing with a problem with N unknowns, the coefficient matrices H and G are dense and non-symmetrical, resulting in O(N ²) arithmetic operations. FMM is used to speed up the solution of the standard boundary element system of equations and reduces the amount of memory required. The core idea behind FMM is to approximate the fundamental solution for BEM using spherical Hankel functions, Legendre polynomials, and plane waves. The coefficient matrices are divided into two portions. The first is the near-field component, which is assessed by integration in the region of the source point. The other is the far-field component, which cannot be calculated directly. Using FMM on a cluster hierarchy decreases the complexity of BEM from O(N ²) to O(N log  N). FMM comes in two varieties. The original FMM (low-frequency technique) is based on the fundamental solution’s series expansion formula, whereas the diagonal form FMM (high-frequency method) is based on the fundamental solution’s plane wave expansion formula. For high-frequency problems, the original FMM is inefficient, while the diagonal form FMM has instability issues when solving low-frequency Helmholtz equations. To circumvent these challenges, the wideband FMM generated by merging the original FMM and the diagonal form FMM can be employed [25]; [26].

The radiated sound power W on an arbitrary closed surface around the radiator may be represented as Equation 12 for radiation into open domains: W A ω = 1 2 ∫ A R p y , ω v f * y , ω d A y , (12)

where

A is a randomly chosen closed surface that encircles the radiator,

p is the sound pressure,

v f * is the conjugate complex of the particle velocity v _f ,

R is the real part of the quantity.

The real component of complex sound power is radiated into the acoustic far field, whereas the imaginary component only contributes to the evanescent near field.

Using BEM to discretize Equation 12 yields a matrix equation for sound power, which is provided by Equation 13: W A ω = 1 2 R p A T S A v A * , (13) S A = ∫ A N f T N f d A y , (14) where

p _A is the nodal sound pressure vector on surface A,

v _A is the particle velocity vector on surface A.

We may readily replace surface A for sound power assessment with the structural surface Γ when the sound power on the structural surface has to be evaluated. As a result, Equation 13 is changed to Equation 15: W ω = 1 2 R p T S v f * . (15)

By resolving Equations 2, 5, 9, respectively, we may obtain the vectors p, u, and v _f in turn. Finally, Equation 15 may be solved to determine the radiated sound power W on the structure surface.

When surface A is not the structural surface, Equation 13 could also be used to determine the sound power W on surface A. The pressure and particle velocity at field point y on surface A in Equation 13 are p _A (y) and v _A (y), respectively. The vectors p _A and v _A may be obtained by solving p _A (y) and v _A (y) at each node on surface A.

Equation 16 may be used to represent the boundary integral equation created on the interaction surface Γ to estimate the sound pressure p _A (y) at a field point y, as follows: p A y = ∫ Γ G x , y q x d Γ x − ∫ Γ F x , y p x d Γ x , (16) G x , y = e i k r 4 π r , (17) where

x is the source point,

y is the field point,

q(x) is the normal derivative of p(x),

F (x, y) is the normal derivative of G (x, y).

Using the continuity condition shown in Equation 18: q A y = − i ω ρ v A y , (18)

Equation 19 is created by differentiating Equation 16 with regard to n(y): v A y = i ω ρ ∂ p A y ∂ n y = i ω ρ ∫ Γ ∂ G x , y ∂ n y q x d Γ x − i ω ρ ∫ Γ ∂ F x , y ∂ n y p x d Γ x . (19)

Discretizing Eqs 16, 19, we get Eqs 20, 21: p A y = g T y q − h T y p , (20) v A y = i ω ρ g 1 T y q − i ω ρ h 1 T y p . (21)

By solving Eqs 20, 21, the nodal sound pressure p _A (y) and the particle velocity v _A (y) are available at each node on surface A. Vectors p _A and v _A on surface A are thus solved. Finally, Equation 13 allows for the solution of the radiated sound power W _A (ω) on surface A.

Equation 8 could be expressed as Equation 22: K − ω 2 M + i ω C − C sf − ω 2 ρ f G S − 1 C fs H u p = B u p = f s p i , (22)

in which we have Equation 23: B = K − ω 2 M + i ω C − C sf − ω 2 ρ f G S − 1 C fs H = A − C sf − ω 2 ρ f G S − 1 C fs H . (23)

By differentiating Equation 22 with regard to the design variable θ, we get Equation 24: ∂ B ∂ θ u p + B ∂ u ∂ θ ∂ p ∂ θ = ∂ f s ∂ θ ∂ p i ∂ θ . (24)

Equations 25, 26 are thus conducted: ∂ B ∂ θ = ∂ A ∂ θ − ∂ C sf ∂ θ − ∂ ω 2 ρ f G S − 1 C fs ∂ θ ∂ H ∂ θ , (25) r = r 1 r 2 = B ∂ u ∂ θ ∂ p ∂ θ = ∂ f s ∂ θ ∂ p i ∂ θ − ∂ B ∂ θ u p . (26)

In the following, several expressions of Eqs 25, 26 are found for different kinds of design variables:

1. When the fluid density ρ _f is chosen to be the design variable θ, Eqs 25, 26 are expressed as Eqs 27, 28:

When dealing with complicated structures, the direct differentiation approach makes it challenging to determine the derivative of A, C _sf , C _fs , S ⁻¹ , H, and G. To overcome this challenge, however, one might use the semi-analytical derivative approach, which uses the finite difference method to compute different coefficient matrices.

Considering the sound power sensitivity on the structure surface Γ, Eq. 33 may be used by differentiating Eq. 15 with regard to design variable θ: ∂ W ∂ θ = 1 2 R ∂ p ∂ θ T w 1 + i ω p T ∂ C fs ∂ θ u * + w 2 T ∂ u ∂ θ * , (33) where

w ₁ = iω C _fs u*,

w 2 T = iω p ^T C _fs .

Introducing the conjugate complex transposed () ^H , we get Eq. 34: R w 2 T ∂ u ∂ θ * = R w 2 H ∂ u ∂ θ . (34)

Applying this to Eq. 33, E. 35 is produced to represent the sound power sensitivity: ∂ W ∂ θ = 1 2 R w 1 T ∂ p ∂ θ + i ω p T ∂ C fs ∂ θ u * + w 2 H ∂ u ∂ θ . (35)

The sum of the first and the third terms in the right side of Eq. 35 can be rewritten as Eq. 36: We may get Eq. 36 by adding the first and third terms of Eq. 35’s right side: w 1 T ∂ p ∂ θ + w 2 H ∂ u ∂ θ = d T ∂ u ∂ θ ∂ p ∂ θ = d T B − 1 r , (36) where d = w 2 * w 1 . (37)

We may represent the sound power sensitivity on the structure surface as Eq. 38 by substituting Eq. 36 into Eq. 35: ∂ W ∂ θ = 1 2 R d T B − 1 r + i ω p T ∂ C fs ∂ θ u * . (38)

Two terms make up the structural surface’s sound power sensitivity. Eq. 38’s first term on the right side can be resolved in one of two ways. One is to first solve B z ̄ = r , a linear system of equations, and then to solve d T z ̄ . The other is to solve zB = d ^T , a linear system of equations, and then zr. Eq. 26 shows that, in contrast to d, r depends on the derivatives of specific terms with respect to the design variable θ. Consequently, m right hand sides r will be produced by m design variables θ _j with j = 1, 2, … , m. The linear system of equations B z ̄ = r must be solved m times using the first method, which takes too much time. The linear system of equations zB = d ^T only has to be solved once for various design variables when using the second method, though.

When an iterative solver, such as GMRES, is used to solve the adjoint equation zB = d ^T , the convergence is low. Eq. 39 may be used to rewrite the adjoint equation: z B = z s z f B = d T , (39) wheres is the degree of freedom of the structure.f is the degree of freedom of the fluid.Equation 40 illustrates a practical way by splitting the adjoint equation into two reduced coupled sensitivity equations: z s A − ω 2 ρ z f G S − 1 C f s = w 2 H , (40) z f H − z s C s f = w 1 T . (41)

Equation 42 is obtained by transforming Equationg 40 into Equationg 41 and removing the vector z _s : z f H − z f G W C sf = w 1 T + w 2 H A − 1 C s f . (42)

The reduced coupled sensitivity equation mentioned above, which is the same as solving Eq. 9, may be used to determine the unidentified fluid vector z _f . The unknown structural vector z _s can then be found using Eq. 40. The second term in Eq. 38 vanishes when the design variable is the fluid density ρ, structural density ρ _s , Poisson’s ratio v, Young’s modulus E, or structural thickness h because ∂ C fs ∂ θ = 0 . However, the term ∂ C fs ∂ θ does not disappear when the structural form parameter is used as the design variable, such as the radius of the spherical shell r. Although difficult, it is feasible to accurately characterize the sensitivity of C _fs analytically. However, using finite differences provides a straightforward and practical solution to this issue.

Equation 43 is created by differentiating Eq. 13 with regard to the design variable θ: ∂ W A ω ∂ θ = 1 2 R ∂ p A ∂ θ T S A v A * + p A T ∂ S A ∂ θ v A * + p A T S A ∂ v A ∂ θ * . (43)

All of the items in vectors ∂ p A ∂ θ and ∂ v A ∂ θ can be solved using the equations of ∂ p A ( y ) ∂ θ and ∂ v A ( y ) ∂ θ . Eqs 44, 45 may be obtained by differentiating Eqs 16, 19: ∂ p A y ∂ θ = ∫ Γ ∂ G x , y ∂ θ q x d Γ x + ∫ Γ G x , y ∂ q x ∂ θ d Γ x + ∫ Γ G x , y q x ∂ d Γ x ∂ θ − ∫ Γ ∂ F x , y ∂ θ p x d Γ x − ∫ Γ F x , y ∂ p x ∂ θ d Γ x − ∫ Γ F x , y p x ∂ d Γ x ∂ θ , (44)

and ∂ v A y ∂ θ = i ω ρ ∫ Γ ∂ 2 G x , y ∂ n y ∂ θ q x d Γ x + i ω ρ ∫ Γ ∂ G x , y ∂ n y ∂ q x ∂ θ d Γ x + i ω ρ ∫ Γ ∂ G x , y ∂ n y q x ∂ d Γ x ∂ θ − i ω ρ ∫ Γ ∂ 2 F x , y ∂ n y ∂ θ p x d Γ x − i ω ρ ∫ Γ ∂ F x , y ∂ n y ∂ p x ∂ θ d Γ x − i ω ρ ∫ Γ ∂ F x , y ∂ n y p x ∂ d Γ x ∂ θ + i ω ∂ ρ − 1 ∂ θ ∫ Γ ∂ G x , y ∂ n y q x d Γ x − i ω ∂ ρ − 1 ∂ θ ∫ Γ ∂ F x , y ∂ n y p x d Γ x , (45)

where ∂ G x , y ∂ θ = − e ikr 4 π r 2 1 − i k r ∂ r ∂ x i ∂ x i ∂ θ − ∂ y i ∂ θ , (46) ∂ F x , y ∂ θ = e ikr 4 π r 3 3 − 3 i k r − k 2 r 2 ∂ r ∂ n x ∂ r ∂ x i − 1 − i k r n i x × ∂ x i ∂ θ − ∂ y i ∂ θ − e ikr 4 π r 2 1 − i k r ∂ r ∂ x i ∂ n i x ∂ θ , (47) ∂ 2 G x , y ∂ n y ∂ θ = e ikr 4 π r 3 3 − 3 i k r − k 2 r 2 ∂ r ∂ n y ∂ r ∂ y i − 1 − i k r n i y × ∂ y i ∂ θ − ∂ x i ∂ θ − e ikr 4 π r 2 1 − i k r ∂ r ∂ y i ∂ n i y ∂ θ , (48) ∂ 2 F x , y ∂ n y ∂ θ = − e ikr 4 π r 4 15 − 15 i k r − 6 k 2 r 2 + i k 3 r 3 ∂ r ∂ n x ∂ r ∂ n y ∂ r ∂ y k − 3 − 3 i k r − k 2 r 2 n k y ∂ r ∂ n x − n k x ∂ r ∂ n y − n j x n j y ∂ r ∂ y k × ∂ y k ∂ θ − ∂ x k ∂ θ + e ikr 4 π r 3 [ 3 − 3 i k r − k 2 r 2 × ∂ r ∂ n y ∂ r ∂ x i ∂ n i x ∂ θ + ∂ r ∂ n x ∂ r ∂ y j ∂ n j y ∂ θ + 1 − i k r ∂ n j x ∂ θ n j y + n j x ∂ n j y ∂ θ ] , (49)

and ∂ d Γ ∂ θ = ∂ 2 x i ∂ θ ∂ x i − ∂ 2 x i ∂ θ ∂ x j n i x n j x d Γ , (50) where ∂ x i ∂ θ will be determined once the design variable has been used to fully parameterize the examined domain’s boundary.

Equations 51 and 52 may be obtained by discretizing Eqs 44, 45: ∂ p A y ∂ θ = g 2 T y , θ q + g T y ∂ q ∂ θ − h 2 T y , θ p − h T y ∂ p ∂ θ , (51) ∂ v A y ∂ θ = i ω ρ g 3 T y , θ q + i ω ρ g 1 T y ∂ q ∂ θ − i ω ρ h 3 T y , θ p − i ω ρ h 1 T y ∂ p ∂ θ + i ω ∂ ρ − 1 ∂ θ g 1 T y q − h 1 T y p , (52) where g ₂, g ₃, h ₂, and h ₃ are coefficient vectors.

Equation 51’s g ₂ and h ₂ as well as Equation 52’s g ₃ and h ₃ disappear when the fluid density ρ is chosen as the design variable. Eqs 51, 52 may be rewritten as Eqs 53, 54 as a result: ∂ p A y ∂ ρ = g T y ∂ q ∂ ρ − h T y ∂ p ∂ ρ , (53) ∂ v A y ∂ ρ = i ω ρ g 1 T y ∂ q ∂ ρ − i ω ρ h 1 T y ∂ p ∂ ρ − i ρ 2 ω g 1 T y q − h 1 T y p . (54)

The variables g ₂ and h ₂ in Equation 51 and g ₃ and h ₃ in Equation 52 disappear when the structural parameter is chosen as the design variable, such as the thickness of the spherical shell, as given in the following numerical example. Eqss 51, 52 are thus equivalent to Eqs 55, 56: ∂ p A y ∂ θ = g T y ∂ q ∂ θ − h T y ∂ p ∂ θ , (55) ∂ v A y ∂ θ = i ω ρ g 1 T y ∂ q ∂ θ − i ω ρ h 1 T y ∂ p ∂ θ . (56) g ₂, g ₃, h ₂, and h ₃ do not disappear when the structural form parameter, such as the radius of the spherical shell, is specified as the design variable. Eqs 51, 52 may be written as Eqs 57, 58 as a result: ∂ p A y ∂ θ = g 2 T y , θ q + g T y ∂ q ∂ θ − h 2 T y , θ p − h T y ∂ p ∂ θ , (57) ∂ v A y ∂ θ = i ω ρ g 3 T y , θ q + i ω ρ g 1 T y ∂ q ∂ θ − i ω ρ h 3 T y , θ p − i ω ρ h 1 T y ∂ p ∂ θ . (58)

The derivatives of p _A (y) and v _A (y) are both shown to be determined by p, q, and their derivatives by Eqs 51, 52. Equations 2, 5, 9, in that order, may be used to produce the vectors p, u, and v _f . Using the continuity condition throughout the interaction surface, the vector q may then be found. We still need to find the solution to the unknown vectors ∂ p ∂ θ and ∂ q ∂ θ , though.

Equation 26 cannot be directly solved because the system matrix is too enormous for problems of this kind. Equations 59 and 60 can be used to separate the system Equation 26: A ∂ u ∂ θ − C sf ∂ p ∂ θ = r 1 , (59) H ∂ p ∂ θ − ω 2 ρ G S − 1 C f s ∂ u ∂ θ = r 2 . (60)

Equation 61 may be obtained by substituting Equation 59 into Equation 60: H ∂ p ∂ θ − G W C sf ∂ p ∂ θ = G W r 1 + r 2 . (61)

Equation 61 and Equation 9 are quite similar, hence the same approach to solving both is used. Equation 59 may be used to solve the unknown vector ∂ u ∂ θ once the foregoing equation has been solved to get the sensitivity of the nodal sound pressure on the structural surface, denoted by the symbol ∂ p ∂ θ . Equation 62 is obtained by differentiating Equation 5 with respect to the design variable and applying the continuity condition over the interaction surface: ∂ q ∂ θ = i ω ∂ ρ ∂ θ v f + ω 2 ρ ∂ S − 1 C f s ∂ θ u + S − 1 C f s ∂ u ∂ θ . (62)

We may find ∂ q ∂ θ by using Equation 62. Equations 51 and 52 may be used to get the derivatives of p _A (y) and v _A (y) at each node on surface A. Equation 63 can be used to describe the derivative of S _A on the right side of Equation 43: ∂ S A ∂ θ = ∫ A N f T N f ∂ d A y ∂ θ . (63)

Once the computing surface does not change when the design variable changes, ∂ S A ∂ θ = 0 causes the second term on the right side of Eq. 43 to disappear. Equation 43 may be used to get the derivative of the radiated sound power on surface A after obtaining the solutions of ∂ p A ∂ θ and ∂ v A ∂ θ .

The sound field of an underwater thin spherical shell centered at position (0, 0, 0) is investigated in this illustration while accounting for a concentrated force F applied at point A (r, 0, 0), where r stands for the radius of the spherical shell, as seen in Figure 1. The material and geometrical features employed in this example are as follows:

radius of the shell is 1.2 m,

thickness of the shell is 0.012 m,

elasticity modulus of the shell is 2.10 × 10¹¹ Pa.

Poisson’s ratio of the shell is 0.3,

structural density is 7.86 × 10³ kg/m³,

fluid density is 1.00 × 10³ kg/m³,

sound velocity in water is 1.482 × 10³ m/s.

The numerical and analytical solutions, expressed in terms of frequencies, for the sound pressure at point (2.4, 0, 0), are shown in Figure 2. The numerical and analytical solutions, expressed in terms of frequencies, for the sound power on the structural shell surface, are shown in Figure 3. The linear systems are solved using the GMRES implementation without preconditioning, and the wideband FMM algorithm is used to speed up the solving process. 6,144 elements make up the discretized thin-shell model. These figures both demonstrate the good agreement between the numerical and analytical solutions, demonstrating that the wideband FMM method preserves the excellent accuracy of traditional BEM.

Figures 4, 5 depicts the structure surface’s sensitivity of sound pressure to the sphere’s radius and its shell thickness, respectively. Figures 6, 7 illustrates the structural surface’s sensitivity of sound power to the sphere’s radius and its shell thickness, respectively. These four figures show a remarkably similar pattern. These figures show that the numerical and analytical results accord rather well. The figures illustrate how the sound pressure or sound power sensitivity is very modest in the low-frequency range but substantially increases at resonance peaks.

The scattering sound field of the underwater submarine model under the influence of plane waves is the subject of this section. The incident wave amplitude is 1.0 Pa, and the plane wave propagates positively along the x-axis. The generic model BeTSSi-Sub presented at the World Digital Simulation Conference in 2002 is adopted by the model. The model has a 0.10 m thickness. The positive x-axis direction is from the bow to the stern, and the origin of the coordinate is located at the intersection of the circle that connects the bow and the hull. For the specific geometric characteristics of the BeTSSi-Sub model, please refer to Figure 3, 4 in [58]. Figure 8 displays the submarine model, with a total of 27,034 elements.

The point for the computation is (40, 0, 0). Figure 9 shows the calculation result of sound pressure at point (40, 0, 0), and Figure 10 is the sensitivity of sound pressure at point (40, 0, 0) to shell thickness. Figures 9, 10 show that when the calculation frequency rises, both the sound pressure and its sensitivity to thickness gradually decline. Given that the sound pressure is significantly higher and more sensitive to changes in structural thickness in the lower frequency band, these two figures show that the lower frequency band is a crucial region for the BeTSSi-Sub model with the current material and geometrical specifications.