An Iterative Physical Acoustics Method for Modeling Acoustic Scattering by Penetrable Objects

1. Introduction

Active acoustic detection is crucial for identifying underwater vehicles, whereby precise modeling of their acoustic scattering characteristics is essential to advancing existing target detection and classification algorithms [1]. Underwater vehicles incorporate diverse structural configurations, encompassing both water-filled structures (e.g., sails and ballast tanks) and pressure-resistant compartments (e.g., equipment bays and crew modules). In high-frequency acoustic analysis of pressure-resistant structures, boundary surfaces are typically treated as wave-impenetrable interfaces through the implementation of hard, soft, or impedance boundary conditions [2]. However, the impenetrable boundary fails to capture the coupled fluid-structure interactions for water-filled thin shells, where acoustic waves not only penetrate the elastic shell but also propagate through the internal fluid, inducing elastic wave excitation [3]. The multiple reflections and transmissions within water-filled shells generate discriminative acoustic scattering features [4], thereby establishing a physical basis for developing feature-driven classification algorithms in underwater target recognition. Therefore, there exists a need to develop computationally efficient methods for characterizing acoustic scattering from water-filled thin shells in order to forego exhaustive experiments and time-consuming simulations.

Many approaches are proposed in the literature to address acoustic scattering problems involving fluid-filled thin shells, and they can be classified into rigorous analytic, exact numerical, and approximate methods. The computation of the interaction between acoustic waves and an elastic shell requires solving a boundary value problem that includes a scalar Helmholtz equation for the acoustic pressure in the fluid medium and a vector Navier equation for the three-dimensional displacements in the elastic medium, which are coupled through boundary conditions at the fluid–elastic interface. The equations can be solved by separation of variables, which can be available to obtain rigorous solutions for scattered fields of simple-shape targets, such as elastic spherical shells [5,6], cylindrical shells [7], and so on. These solutions for canonical shapes, despite their limited practicality, establish critical foundations for computational validation or mechanism analysis in acoustic scattering. The governing equations in the elastic medium can be formulated either in differential form, which is typically approximated using the Finite Element Method (FEM) [8], or in integral form, which is approximated using the Boundary Element Method (BEM) [9]. The coupling of the interior problem for the elastic shell with the exterior problem in the fluid leads to the development of hybrid numerical approaches, specifically the FEM/BEM [10,11]. Although they can provide numerically exact results, their computational efficiency significantly diminishes at higher frequencies due to increasing computation time and storage demands.

To improve computational efficiency in modeling acoustic scattering from elastic shell structures, it is beneficial to replace direct numerical simulations of elastoacoustic interactions with simplified theoretical approaches based on specific assumptions. The main assumptions are: 1) the thickness of the shell is much smaller than the linear size of the structure, and 2) the thickness is much smaller than the acoustic wavelength in the fluid. Under these conditions, it is logical to apply established theories for elastic plates and shells [12] as a practical framework. A shell can be approximated by introducing effective boundary conditions, effectively reducing the problem from a vector-based formulation to a simpler scalar framework. In the simplest case, this condition can express Newton's second law, and the shell can be considered purely inertial [13,14]. Further, to account for the complex nature of the interaction of acoustic waves with elastic thin shells, some kind of effective boundary conditions are derived in the form of surface differential equations [15]. In the end, the problem is reduced to a system of hypersingular integral equations, and then solved by the BEM. However, in the middle and higher frequency range, the computing resources for larger thin shell structures are usually not available in practice.

For modeling mid-to-high frequency acoustic scattering from underwater complex shell structures, approximate methodologies such as ray-based approaches and physical acoustics (PA)-based methods are helpful alternatives to full-wave numerical simulations. The ray synthesis [16,17] decomposes the scattering phenomenon into geometric and elastic wave components, which are subsequently synthesized via coherent superposition. While offering a clear physical interpretation of wave interaction mechanisms, its application remains limited to simple geometries due to challenges in calculating elastic scattering from complex shell structures [18]. The geometric ray acoustics [19], applied in practical engineering analyses of complex shell structures, enables rapid identification of echo highlights through superposition of specular reflection components. However, it fails to account for the contributions of non-specular and elastic scattering. PA-based methods are applicable for analyzing both impenetrable boundaries [20,21] and penetrable interfaces [22,23,24] in acoustic wave propagation. A recent approach integrates geometric ray theory with physical acoustics to predict the target strength (TS) of water-filled thin-shell structures [25]. By incorporating PA, this hybrid approach accounts for non-specular geometric echoes, which are essential for characterizing complex scattering mechanisms. However, the approach limits its analysis to the first-order reflection and transmission, neglecting higher-order reflections and transmissions, which dominate in penetrable shells.

This study aims to develop a simple, computationally efficient, and physically intuitive approximate model for analyzing the acoustic scattering characteristics of arbitrarily shaped water-filled thin shells immersed in a medium. Building upon well-established theories of elastic thin plates [26], this work extends the original iterative physical acoustics (IPA) method [27,28] to the case of penetrable thin-shell structures. The proposed method relies on two fundamental assumptions: (i) the thin-shell hypothesis, which ensures inertial dominance over elastic effects across a broad frequency range [14], and (ii) the Kirchhoff approximation, which governs reflection/transmission behavior via localized plane-wave principles [29]. The method determines the complex reflection/transmission coefficients at the incident surface using local effective boundary conditions, initializes surface pressure distributions via scalar formalism, and iteratively refines these fields to convergence. This iterative procedure inherently captures the interplay of frequency-dependent multiple reflections and transmissions through the structure. Leveraging a prior discretization framework [27,28,30], an extended IPA-based facet scattering formulation for thin shells can be derived, balancing computational efficiency with physical fidelity to capture multiple reflections and transmissions.

The remainder of this work is organized as follows. Section 2 briefly introduces the Helmholtz integral equation, Kirchhoff approximation for reflection and transmission, and their implementation in modeling acoustic scattering from thin-shell structures. Then, the discretization scheme and numerical implementation of the thin-shell integral formulation are presented. Numerical validation and simulations are described using representative cases. In Section 3, the proposed method is applied to model the acoustic multiple reflection and transmission from a scaled sail of Benchmark Target Strength Simulations (BeTSSi IIB) [31], and an acoustic scattering experiment on the scaled sail is conducted and analyzed. Finally, Section 4 presents the concluding remarks.

2. Methodology

2.1. Thin-shell integral formulation

Consider a thin shell submerged in an infinite acoustic medium with a speed of sound $c_0$ and mean density ${\rho }_0$, as shown in Figure 1. The shell is made of an elastic material characterized by a mass density ${\rho }_e$, Young's modulus $E$, and Poisson's ratio $\nu$, and has a thickness $h$. The neutral surface of the thin shell is denoted by $S$, and an imaginary surface $s$ is constructed to divide the acoustic domain. The shell divides two subdomains, denoted by the exterior subdomain ${\Omega }^{+}$, bounded by the outer boundary $S^{+}$, and the interior subdomain ${\Omega }^{-}$, bounded by the inner boundary $S^{-}$. $n$ is the normal vector pointing in the direction of ${\Omega }^{+}$. The shell is exposed to an incident wave pressure $p^{inc}$ in a scattering problem. The surface pressure on the $S^{+}$ is denoted by $p^{+}$, and that on the $S^{-}$ is denoted by $p^{-}$. By applying the Helmholtz integral equation to each subdomain, the thin-shell integral equation [32] can be derived and expressed in the following form

$$\begin{array}{l}{p}^{+}\left(r\right)+p^{-}\left(r\right)=2p^{inc}\left(r\right)+2{\displaystyle \underset{S}{\int }\left\{\left[p^{+}\left(r_s^{\text{'}}\right)-p^{-}\left(r_s^{\text{'}}\right)\right]\right.}{\displaystyle \frac{\partial G\left(\left.r\right|r_s^{\text{'}}\right)}{\partial n\left(r_s^{\text{'}}\right)}}\\ -\left[{\displaystyle \frac{\partial p^{+}\left(r_s^{\text{'}}\right)}{\partial n\left(r_s^{\text{'}}\right)}}-{\displaystyle \frac{\partial p^{-}\left(r_s^{\text{'}}\right)}{\partial n\left(r_s^{\text{'}}\right)}}\right]\left.G\left(\left.r\right|r_s^{\text{'}}\right)\right\}dS\text{'},r\in S,\end{array}$$

(1)

$$\begin{array}{l}{p}^{scat}\left(r\right)={\displaystyle \underset{S}{\int }\left\{\left[p^{+}\left(r_s\right)-p^{-}\left(r_s\right)\right]\right.}{\displaystyle \frac{\partial G\left(r\left|r_s\right.\right)}{\partial n\left(r_s\right)}}\\ -\left[{\displaystyle \frac{\partial p^{+}\left(r_s\right)}{\partial n\left(r_s\right)}}-{\displaystyle \frac{\partial p^{-}\left(r_s\right)}{\partial n\left(r_s\right)}}\right]\left.G\left(r\left|r_s\right.\right)\right\}dS,r\notin S.\end{array}$$

(2)

where $r_s$ and $r_s^{\text{'}}$ are the distinct position vectors of arbitrary points on $S$; $p^{scat}$ is the scattered pressure at an arbitrary spatial coordinate $r$ outside $S$; $G\left(r\left|r_s\right.\right)=e^{ik\left|r-r_s\right|}/4\pi \left|r-r_s\right|$ is the free-field Green's function, and $k$ is the medium wavenumber.

Eq. (1) itself is not sufficient for solving the problem ascribed to the additional unknown functions in Eq. (2). One more equation is then required to supplement Eq. (1) for the solution. The normal derivative of Eq. (1) can be written as

(3)

The iterative solutions of Eqs. (1) and (3) can be represented in the form

$$\begin{array}{c}{\mu }_q^{+}\left(r_s\right)+{\mu }_q^{-}\left(r_s\right)=2p^{inc}\left(r\right)+2{\displaystyle \underset{S}{\int }\left\{\left[{\mu }_{q-1}^{+}\left(r_s^{\text{'}}\right)-{\mu }_{q-1}^{-}\left(r_s^{\text{'}}\right)\right]\right.}{\displaystyle \frac{\partial G\left(\left.r_s\right|r_s^{\text{'}}\right)}{\partial n\left(r_s^{\text{'}}\right)}}\\ -\left[{\displaystyle \frac{\partial {\mu }_{q-1}^{+}\left(r_s^{\text{'}}\right)}{\partial n\left(r_s^{\text{'}}\right)}}-{\displaystyle \frac{\partial {\mu }_{q-1}^{-}\left(r_s^{\text{'}}\right)}{\partial n\left(r_s^{\text{'}}\right)}}\right]\left.G\left(\left.r_s\right|r_s^{\text{'}}\right)\right\}dS\text{'},\end{array}$$

(4a)

$${\mu }_0^{\pm }\left(r_s\right)=0,p^{+}\left(r_s\right)=\underset{q\to \infty }{\lim }{\mu }_q^{+}\left(r_s\right),p^{-}\left(r_s\right)=\underset{q\to \infty }{\lim }{\mu }_q^{-}\left(r_s\right),$$

(4b)

$$\begin{array}{c}{\displaystyle \frac{\partial {\mu }_q^{+}\left(r_s\right)}{\partial n\left(r_s\right)}}+{\displaystyle \frac{\partial {\mu }_q^{-}\left(r_s\right)}{\partial n\left(r_s\right)}}=2{\displaystyle \frac{\partial p^{inc}\left(r\right)}{\partial n\left(r_s\right)}}+2{\displaystyle \underset{S}{\int }\left\{\left[{\mu }_{q-1}^{+}\left(r_s^{\text{'}}\right)-{\mu }_{q-1}^{-}\left(r_s^{\text{'}}\right)\right]\frac{{\partial }^{2}G\left(r_s\left|r_s^{\text{'}}\right.\right)}{\partial n\left(r_s\right)\partial n\left(r_s^{\text{'}}\right)}\right.}\\ -\left[{\displaystyle \frac{\partial {\mu }_{q-1}^{+}\left(r_s^{\text{'}}\right)}{\partial n\left(r_s^{\text{'}}\right)}}-{\displaystyle \frac{\partial {\mu }_{q-1}^{-}\left(r_s^{\text{'}}\right)}{\partial n\left(r_s^{\text{'}}\right)}}\right]\left.{\displaystyle \frac{\partial G\left(r_s\left|r_s^{\text{'}}\right.\right)}{\partial n\left(r_s\right)}}\right\}dS\text{'},\end{array}$$

(5a)

$${\displaystyle \frac{\partial {\mu }_0^{\pm }\left(r_s\right)}{\partial n\left(r_s\right)}}=0,{\displaystyle \frac{\partial p^{+}\left(r_s\right)}{\partial n\left(r_s\right)}}=\underset{q\to \infty }{\lim }{\displaystyle \frac{\partial {\mu }_q^{+}\left(r_s\right)}{\partial n\left(r_s\right)}},{\displaystyle \frac{\partial p^{-}\left(r_s\right)}{\partial n\left(r_s\right)}}=\underset{q\to \infty }{\lim }{\displaystyle \frac{\partial {\mu }_q^{-}\left(r_s\right)}{\partial n\left(r_s\right)}}.$$

(5b)

Thus, in the $q$th approximation, the surface pressure and its normal derivative are generated by the incident wave, surface pressure, and its normal derivative originating as a result of the $\left(q-1\right)$th iteration.

Denoting the plane-wave reflection and transmission coefficients as ${\Re }_s$ and ${\mathfrak{T}}_s$, the initial surface acoustic field ($q=1$) and its normal derivative on the reflection surface $S^{+}$ are given in the single scattering Kirchhoff approximation as [33]

$${\mu }_1^{+}\left(r_s\right)=\left(1+{\Re }_s\right)p^{inc},$$

(6a)

$${\displaystyle \frac{\partial {\mu }_1^{+}\left(r_s\right)}{\partial n\left(r_s\right)}}=\left(1-{\Re }_s\right){\displaystyle \frac{\partial p^{inc}}{\partial n\left(r_s\right)}},$$

(6b)

and the initial surface acoustic field ($q=1$) and its normal derivative on the transmission surface $S^{-}$ can be expressed as [23]

$${\mu }_1^{-}\left(r_s\right)={\mathfrak{T}}_s{p}^{inc},$$

(7a)

$${\displaystyle \frac{\partial {\mu }_1^{-}\left(r_s\right)}{\partial n\left(r_s\right)}}={\mathfrak{T}}_s{\displaystyle \frac{\partial p^{inc}}{\partial n\left(r_s\right)}}.$$

(7b)

The relationship

$${\mu }_1^{+}\left(r_s\right)-{\mu }_1^{-}\left(r_s\right)=\left(1+{\Re }_s-{\mathfrak{T}}_s\right)p^{inc},$$

(8a)

$${\displaystyle \frac{\partial {\mu }_1^{+}\left(r_s\right)}{\partial n\left(r_s\right)}}-{\displaystyle \frac{\partial {\mu }_1^{-}\left(r_s\right)}{\partial n\left(r_s\right)}}=\left(1-{\Re }_s-{\mathfrak{T}}_s\right){\displaystyle \frac{\partial p^{inc}}{\partial n\left(r_s\right)}},$$

(8b)

motivates us to rewrite Eqs. (4) and (5) through the application of the Kirchhoff approximation at each iteration as

$$\begin{array}{l}{\mu }_q^{+}\left(r_s\right)-{\mu }_q^{-}\left(r_s\right)=\left(1+{\Re }_s-{\mathfrak{T}}_s\right){\displaystyle \underset{S}{\int }\left\{\frac{\partial G\left(\left.r_s\right|r_s^{\text{'}}\right)}{\partial n\left(r_s^{\text{'}}\right)}\right.}\left[{\mu }_{q-1}^{+}\left(r_s^{\text{'}}\right)\right.\left.-{\mu }_{q-1}^{-}\left(r_s^{\text{'}}\right)\right]\\ -G\left(\left.r_s\right|r_s^{\text{'}}\right)\left.\left[{\displaystyle \frac{\partial {\mu }_{q-1}^{+}\left(r_s^{\text{'}}\right)}{\partial n\left(r_s^{\text{'}}\right)}}-{\displaystyle \frac{\partial {\mu }_{q-1}^{-}\left(r_s^{\text{'}}\right)}{\partial n\left(r_s^{\text{'}}\right)}}\right]\right\}dS\text{'},q\ge 2\end{array}$$

(8c)

$$\begin{array}{l}{\displaystyle \frac{\partial {\mu }_q^{+}\left(r_s\right)}{\partial n\left(r_s\right)}}-{\displaystyle \frac{\partial {\mu }_q^{-}\left(r_s\right)}{\partial n\left(r_s\right)}}=\left(1-{\Re }_s-{\mathfrak{T}}_s\right){\displaystyle \underset{S}{\int }\left\{\frac{{\partial }^{2}G\left(r_s\left|r_s^{\text{'}}\right.\right)}{\partial n\left(r_s\right)\partial n\left(r_s^{\text{'}}\right)}\right.}\left[{\mu }_{q-1}^{+}\left(r_s^{\text{'}}\right)\right.\left.-{\mu }_{q-1}^{-}\left(r_s^{\text{'}}\right)\right]\\ \left.-{\displaystyle \frac{\partial G\left(r_s\left|r_s^{\text{'}}\right.\right)}{\partial n\left(r_s\right)}}\left[{\displaystyle \frac{\partial {\mu }_{q-1}^{+}\left(r_s^{\text{'}}\right)}{\partial n\left(r_s^{\text{'}}\right)}}-{\displaystyle \frac{\partial {\mu }_{q-1}^{-}\left(r_s^{\text{'}}\right)}{\partial n\left(r_s^{\text{'}}\right)}}\right]\right\}dS\text{'},q\ge 2\end{array}$$

(8d)

Substituting Eq. (8) into Eq. (2) yields

$$\begin{array}{c}{p}^{scat}\left(r\right)={\displaystyle \underset{S}{\int }\left\{{\displaystyle {\sum }_{q=1}^Q\left[{\mu }_q^{+}\left(r_s\right)-{\mu }_q^{-}\left(r_s\right)\right]}\right.}{\displaystyle \frac{\partial G\left(r_B\left|r_s\right.\right)}{\partial n\left(r_s\right)}}\\ -{\displaystyle {\sum }_{q=1}^Q\left[\frac{\partial {\mu }_q^{+}\left(r_s\right)}{\partial n\left(r_s\right)}-\frac{\partial {\mu }_q^{-}\left(r_s\right)}{\partial n\left(r_s\right)}\right]}\left.G\left(r_B\left|r_s\right.\right)\right\}dS,\end{array}$$

(9)

where $Q$ denotes the number of iterations. The convergence of Eq. (9) can be effectively monitored using a suitable global parameter given by

$$\left\{\begin{array}{l}{‖\left[{\mu }_q^{+}\left(r_s\right)-{\mu }_q^{-}\left(r_s\right)\right]-\left[{\mu }_{q-1}^{+}\left(r_s\right)-{\mu }_{q-1}^{-}\left(r_s\right)\right]‖}_2/{‖{\mu }_1^{+}\left(r_s\right)-{\mu }_1^{-}\left(r_s\right)‖}_2<\sigma ,\\ {‖\left[{\displaystyle \frac{\partial {\mu }_q^{+}\left(r_s\right)}{\partial n\left(r_s\right)}}-{\displaystyle \frac{\partial {\mu }_q^{-}\left(r_s\right)}{\partial n\left(r_s\right)}}\right]-\left[{\displaystyle \frac{\partial {\mu }_{q-1}^{+}\left(r_s\right)}{\partial n\left(r_s\right)}}-{\displaystyle \frac{\partial {\mu }_{q-1}^{-}\left(r_s\right)}{\partial n\left(r_s\right)}}\right]‖}_2/{‖{\displaystyle \frac{\partial {\mu }_1^{+}\left(r_s\right)}{\partial n\left(r_s\right)}}-{\displaystyle \frac{\partial {\mu }_1^{-}\left(r_s\right)}{\partial n\left(r_s\right)}}‖}_2<\sigma .\end{array}\right.$$

(10)

From Eqs. (8) and (9), the scattered pressure field is fundamentally governed by the interfacial discontinuities in both surface pressure and its normal gradient across the shell boundary. These governing physical quantities, in turn, are determined by the inherent properties of the shell. Assuming that the Helmholtz numbers $ka\gg 1$ and $h/a\ll 0.05$ ($a$ being the characteristic size of the scattering object), the quantities ${\Re }_s$ and ${\mathfrak{T}}_s$ can be derived by the approximate theory of vibration of an infinite elastic thin plate as [26]

$$\Re ={\displaystyle \frac{Z_p^M}{Z_p^M+2Z_p^a}},\mathfrak{T}={\displaystyle \frac{2Z^a}{Z_p^M+2Z_p^a}},$$

(11)

where $Z_p^M$ and $Z_p^a$ are the mechanical impedance of a thin plate and the acoustic impedance of fluid, respectively, as

$$Z_p^M={\displaystyle \frac{iD}{\omega }}\left[{\left(k\sin {\theta }_i\right)}^4-{\gamma }^4\right],Z_p^a={\rho }_{0}c/\cos {\theta }_i,$$

(12)

where $D=E{h}^3/12\left(1-{\nu }^2\right)$ and ${\gamma }^4={\omega }^2{\rho }_{e}h/D$, and ${\theta }_i$ is the incidence angle. Substitution of Eq. (11) into Eqs. (8) and (9) yields a simplified formulation as

$$p^{scat}\left(r\right)={\displaystyle \underset{S}{\int }{\displaystyle {\sum }_{q=1}^Q\left[{\mu }_q^{+}\left(r_s\right)-{\mu }_q^{-}\left(r_s\right)\right]}}{\displaystyle \frac{\partial G\left(r_B\left|r_s\right.\right)}{\partial n\left(r_s\right)}}dS,$$

(13a)

$$\left\{\begin{array}{l}{\mu }_1^{+}\left(r_s\right)-{\mu }_1^{-}\left(r_s\right)=\left(1+{\Re }_s-{\mathfrak{T}}_s\right)p^{inc},q=1\\ {\mu }_q^{+}\left(r_s\right)-{\mu }_q^{-}\left(r_s\right)=\left(1+{\Re }_s-{\mathfrak{T}}_s\right){\displaystyle \underset{S}{\int }\frac{\partial G\left(\left.r_s\right|r_s^{\text{'}}\right)}{\partial n\left(r_s^{\text{'}}\right)}}\left[{\mu }_{q-1}^{+}\left(r_s^{\text{'}}\right)\right.\left.-{\mu }_{q-1}^{-}\left(r_s^{\text{'}}\right)\right]dS\text{'},q\ge 2\end{array}\right.$$

(13b)

As the plate thickness approaches the thin-plate limit, the integral component of the single-layer potential asymptotically diminishes. For the rigid target scattering, i.e., ${\Re }_s=1$ and ${\mathfrak{T}}_s=0$, Eq. (13) reduces to

$$p^{scat}\left(r\right)={\displaystyle \underset{S}{\int }{\displaystyle {\sum }_{q=1}^Q{\mu }_q\left(r_s\right)}}{\displaystyle \frac{\partial G\left(r_B\left|r_s\right.\right)}{\partial n\left(r_s\right)}}dS,$$

(14a)

$$\left\{\begin{array}{l}{\mu }_1\left(r_s\right)=2p^{inc},q=1\\ {\mu }_q\left(r_s\right)=2{\displaystyle \underset{S}{\int }{\mu }_{q-1}\frac{\partial G\left(\left.r_s\right|r_s^{\text{'}}\right)}{\partial n\left(r_s^{\text{'}}\right)}}dS\text{'},q\ge 2\end{array}\right.$$

(14b)

The integral formulation for the scattered pressure field generated by a thin shell structure at arbitrary spatial coordinates, which accounts for multiple reflected and transmitted wave components, is established in Eq. (13). The expression specifically characterizing the multiply reflected waves of a rigid target, as derived in Eq. (14), has been validated in previous studies [27,28].

2.2. Discretization and numerical implementation

A numerical solution of Eq. (13) can be achieved by discretizing the thin-shell neutral surface S into $N$ triangular facets. The scattered pressure of the thin shell at the field point $B$ can be approximated as the coherent superposition of $Q$ orders of reflected and transmitted fields generated by $N$ facets. Figure 2 is a sketch of two facets $m$ and $n$ of the problem. The notation used is as follows:

1) $r_A$ is the source position vector from an origin point $O$ to a source point $A$, $r_B$ is the field position vector from $O$ to the field point $B$;

2) $r_n$ is the position vector from $O$ to an arbitrary point $F_n$ on the $n$th facet's surface $S_n$, $r_{\mathrm{c}n}$ is the position vector from $O$ to a centroid $C_n$ on $S_n$, and ${\sigma }_n=r_n-r_{\mathrm{c}n}$;

3) $R_{An}=r_A-r_{\mathrm{c}n}$, $R_{An}=\left|R_{An}\right|$, and ${\widehat{R}}_{An}=R_{An}/R_{An}$; $R_{Bn}=r_B-r_{\mathrm{c}n}$, $R_{Bn}=\left|R_{Bn}\right|$, and ${\widehat{R}}_{Bn}=R_{Bn}/R_{Bn}$;

4) ${\widehat{n}}_n=\widehat{n}\left(r_{\mathrm{c}n}\right)$ is an outward unit normal to $S_n$ at the point $C_n$;

5) Note that the subscript "$m$" can substitute for "$n$" to denote the above definitions regarding the $m$th facet. Also, $R_{nm}=r_{\mathrm{c}n}-r_{\mathrm{c}m}$, $R_{nm}=\left|R_{nm}\right|$, and ${\widehat{R}}_{nm}=R_{nm}/R_{nm}$.

Under the discrete framework, Eq. (13) can be recast in a facet-based representation as follows:

$$p^{scat}\left(\left.R_{An}\right|R_{Bn}\right)={\displaystyle \sum _{n=1}^N\left({\displaystyle \sum _{q=1}^Q\delta {\mu }_{qn}}\right)\partial G_n},$$

(15a)

$$\left\{\begin{array}{l}q=1,\delta {\mu }_{1n}=\left(1+{\Re }_n-{\mathfrak{T}}_n\right){\xi }_{An}{\displaystyle \frac{e^{ik{R}_{An}}}{R_{An}}},\\ q\ge 2,\delta {\mu }_{qn}=\left(1+{\Re }_n-{\mathfrak{T}}_n\right){\displaystyle \sum _{m=1\left(\ne n\right)}^N\delta {\mu }_{\left(q-1\right)m}}\partial G_{nm}.\end{array}\right.$$

(15b)

with

$$\partial G_n=-{\widehat{n}}_n\cdot {\widehat{R}}_{Bn}{I}_{Bn}{e}^{ik{R}_{Bn}}{\displaystyle \frac{ik{R}_{Bn}-1}{{R_{Bn}}^2}},$$

(15c)

$$\partial G_{nm}=\left\{\begin{array}{l}-{\zeta }_{nm}{\widehat{n}}_m\cdot {\widehat{R}}_{nm}{I}_{nm}{e}^{ik{R}_{nm}}{\displaystyle \frac{ik{R}_{nm}-1}{{R_{nm}}^2}},m\ne n\\ 0,m=n,\end{array}\right.$$

(15d)

where $\delta {\mu }_{qn}={\mu }_q^{+}\left(r_n\right)-{\mu }_q^{-}\left(r_n\right)$ is the $q$th order of jump of surface pressure at the point $C_n$ on $S_n$, and $\delta {\mu }_{\left(q-1\right)m}={\mu }_{q-1}^{+}\left(r_m\right)-{\mu }_{q-1}^{-}\left(r_m\right)$ is the $\left(q-1\right)$th order of jump of surface pressure at the point $C_m$ on $S_m$, respectively; ${\xi }_{An}$ is the incident field visibility function and represents the deterministic quantities that the centroid at $r_{\mathrm{c}n}$ is visible to the source at $r_A$; ${\zeta }_{nm}$ is the mutual visibility function and represents the quantities that one centroid at $r_{\mathrm{c}n}$ is visible to another centroid at $r_{\mathrm{c}m}$; $I_{Bn}$ and $I_{nm}$ are surface integral kernels in terms of the $n$th and $m$th facets, respectively, which can be evaluated using Gordon's integral algorithm [34]. These auxiliary functions are as follows:

$${\xi }_{An}=\xi \left(r_A\left|r_{cn}\right.\right)=\left\{\begin{array}{l}1,{\widehat{R}}_{An}\cdot {\widehat{n}}_n>0,\\ 0,\mathrm{otherwise},\end{array}\right.$$

(16a)

$${\zeta }_{nm}=\zeta \left(r_{cm}\left|r_{cn}\right.\right)=\left\{\begin{array}{l}1,{\widehat{R}}_{nm}\cdot {\widehat{n}}_m>0,{\widehat{R}}_{nm}\cdot {\widehat{n}}_n<0,\\ 0,\mathrm{otherwise},\end{array}\right.$$

(16b)

$$I_{Bn}=I\left({\widehat{R}}_{Bn}\left|{\sigma }_n\right.\right)={\displaystyle {\int }_{S_n}{e}^{-ik{\widehat{R}}_{Bn}\cdot {\sigma }_n}}dS,$$

(17a)

$$I_{nm}=I\left({\widehat{R}}_{nm}\left|{\sigma }_m\right.\right)={\displaystyle {\int }_{S_m}{e}^{-ik{\widehat{R}}_{nm}\cdot {\sigma }_m}}dS.$$

(17b)

Further details regarding the derivation of Eq. (15) are available in Refs. [27,28], which provide a comprehensive treatment of the discretization methodology and numerical implementation of Eq. (14).

2.3. Numerical validation

This section focuses on the canonical problem of acoustic scattering from a water-filled spherical shell immersed in an infinite fluid domain, as illustrated in Figure 3. We consider a plane wave incident on the spherical shell in the positive $X$ direction, with all geometric and material parameters provided in Table 1. This well-established benchmark case, for which exact analytical solutions are available [5,6], provides rigorous validation of the IPA and its numerical implementation. Unless otherwise stated, a triangular mesh size of $\lambda /6$ and a convergence threshold of $\sigma ={10}^{-3}$ are set for the IPA simulations in this work. All simulations are performed on a personal computer with the Intel Xeon Pentium 8375C CPU with 3.5 GHz and 1.0 TB RAM.

A comparison between IPA results and analytical solutions is conducted across varying Helmholtz numbers $ka$ and shell thicknesses $h$. Monostatic target strength (TS)[2] results for $h=1$, 3, 10 and 20 mm are respectively presented in Figure 4a–d. As theoretically anticipated, the prediction accuracy of TS increases with $ka$. At elevated $ka$regimes, the predicted TS closely matches analytical solutions, while in the low-$ka$ domain, accuracy progressively decreases. Crucially, reduced thickness-to-radius ratio ($h/a$) excites higher-density elastic resonance spectra, yet the characteristic interference patterns remain accurately predicted. Maximum discrepancies emerge in the low-$ka$, low-$h/a$ regime, attributable to unmodeled elastic wave coupling mechanisms. The IPA solutions reproduce the general behavior of the analytical solutions, including the interference of multiple reflected and transmitted wave components, phenomena notably absent in both rigid and pressure-release boundary approximations.

For a water-filled spherical shell with thickness $h=3$ mm, Figure 5a–d presents near-field scattered pressure level distributions at 5 kHz ($ka\approx 20.9$) and 10 kHz ($ka\approx 41.9$), with the upper panels displaying the IPA results and the bottom panels showing the analytical solutions. The computational domain comprises 160,801 field points distributed across a 4.0 m × 4.0 m square region. Figure 6a–b compares the far-field angular distributions of bistatic TS between IPA and analytical solutions at identical frequencies, with the bistatic angle ${\theta }_b$ sampled at 0.5° resolution. Overall, both methods exhibit strong agreement in resolving interference structures and pressure field extrema localization at high $ka$ values. The angular density of interference fringes scales with excitation frequency. Significantly, the proposed methodology provides comprehensive acoustic field prediction capabilities, encompassing internal cavity fields, external near-fields, and far-field scattering patterns within a unified computational framework. Localized discrepancies stem from IPA's inability to model shear wave propagation in the elastic solid.

3. Sail scattering modeling and experimental validation

3.1. IPA-based sail scattering modeling

The BeTSSi IIB model has emerged as a widely utilized benchmark in submarine acoustics research, distinguished by its characteristic thin-shell sail structure that significantly contributes to its scattering signature analysis. Accordingly, a 1:10 geometric scale model of the sail was developed to facilitate comprehensive modeling and experimental validation. As shown in Figure 7, the model incorporates symmetrical airfoil profiles on both its upper and lower boundaries, mathematically defined by the following equation

$$\left\{\begin{array}{l}y=0.5t\left(a_0{x}^{0.5}-a_{1}x-a_2{x}^2+a_3{x}^3-a_4{x}^4\right),0\le x\le 1,\\ a_0=0.2969,a_1=0.1265,a_2=0.3498,a_3=0.2839,a_4=0.1045\end{array}\right.$$

(18)

where $t=0.1626$ at the upper boundary and $t=0.1816$ at the lower boundary, respectively. The principal geometric parameters are summarized in Table 2. The model consists of a water-filled stainless steel shell, with the material parameters as specified in Table 1.

To compare the accuracy and computational efficiency between the FEM and IPA for general thin-shell structures, the sail model in Figure 7 is analyzed. The FEM simulation is implemented using the Acoustic-Solid Interaction, Frequency Domain multiphysics module within COMSOL Multiphysics commercial software [35]. Variations in the monostatic TS against frequency at the incident azimuth angle ${\theta }_i$ of ${90}^{°}$ are illustrated in Figure 8a. The error between the IPA result and the FEM result decreases with increasing frequency. The resonance peaks existing on the FEM curves do not change the overall tendency of the variation. At high frequencies, the TS predicted by the IPA is very close to the FEM. At lower frequencies, the accuracy decreases but the trend remains physically consistent. The corresponding CPU time used for computing TS in the FEM and IPA is given in Figure 8b. It is clear from the plot that the IPA can greatly shorten the computation time compared with the time spent by FEM, with particularly substantial efficiency gains observed at higher frequencies. One should be noted that the FEM exhibits substantial memory requirements, with the computational capacity being constrained for solving higher-frequency problems.

3.2. Sail scattering experiment and validation

To verify the effectiveness of the proposed IPA method and demonstrate scattering characteristics of a thin-shell target, we perform an acoustic scattering experiment on the scaled sail model immersed in an anechoic tank of 30m$\times$15m$\times$10m. To simulate an approximate free-field environment within the tank while ensuring noise isolation, the interiors are lined with sound-absorbing rubber panels designed for optimal sound absorption. As shown in Figure 9a, the experimental stainless-steel model is constructed with the same geometric dimensions and acoustic parameters as those used in the numerical simulation. Figure 9b,c shows the experimental setup. The experimental model spins around its axis at ${0.5}^{°}/\mathrm{s}$ in sequence and is controlled by a mechanical turntable, whereas the transducer and receiver are fixed. The transducer has a frequency range of 60–120 kHz with an average transmitting voltage response of 146 dB re $\mathsf{\mu }\mathrm{Pa}/\mathrm{V}$ re 1 m. The receiver is a piezoelectric hydrophone (Type 8103, Brüel & Kjaer, Nairom, Denmark) operating over a frequency range of 0.1 Hz–180 kHz with a receiving sensitivity of −211 dB re $1\mathrm{V}/\mathsf{\mu }\mathrm{Pa}$. An uninterrupted power supply powers the entire measurement system to remove AC interference.

The transmitted signal is a linear frequency-modulated chirp spanning from 60 kHz to 80 kHz, with a pulse width of 0.2 ms and a burst interval of 500 ms to ensure steady-state measurement conditions. A matched filter is implemented to enhance signal-to-noise ratio and time-delay resolution. Figure 10a,b compares the time-angle spectra of matched filtered echoes obtained from the scanned measurement and IPA simulation, showing close agreement. Both clearly exhibit primary reflections centered at 12.6 ms (outer surface signature) and higher-order reflections and transmissions beyond 12.8 ms across ${15}^{°}$-${165}^{°}$ and ${195}^{°}$-${345}^{°}$ incidences, though experimental features were less distinct due to noise. Minor inconsistencies in echo envelope morphology suggest possible misalignment between the transducer, receiver, and model centroid. As depicted in Figure 10c, the incident wave illuminates the thin shell structure, generating reflected (superscript $"R"$) and transmitted (superscript $"T"$) wave components. The primary reflected wave ($1^R$) propagates away from the lower surface, while the field transmitted into the interior fluid along the upward propagation direction ($1^T$) is reflected onto the upper surface along the downward propagation direction ($2^R$), and then reflected onto the lower surface and so on. This process continues with successive reflections and transmissions at both upper and lower surfaces. The fundamental propagation paths ① to ⑥ in Figure 10c correspond directly to the temporal features highlighted in Figure 10a,b. The high-order highlights associated with longer paths move toward the farther range than the primary reflections. At paths longer than ⑥, the echo signal is too weak to be shown.

Figure 11 compares the monostatic TS of the sail model obtained using the IPA and experimental measurements at 70 kHz and 80 kHz. Experimental and simulated results exhibit great agreement around the main lobes of scattering patterns, demonstrating that the IPA method effectively captures dominant scattering mechanisms and reliably predicts high-frequency scattering fields of penetrable thin shells. However, discrepancies between simulated and experimental TS are evident in the bow and stern angular regions, where simulations underestimate measured values. This deviation primarily arises from two unmodeled physical mechanisms: (i) The increased curvature of the bow shell triggers the excitation of multimodal surface elastic waves—particularly Franz waves—while enhancing radiative efficiency of Rayleigh and Stoneley waves [36], collectively amplifying scattered energy in experimental measurements; with a secondary contribution from (ii) edge diffraction echoes induced by the sharp trailing-edge discontinuity at the stern [37].

4. Conclusions

This study proposes an IPA method for efficient modeling of acoustic scattering from thin shells immersed in fluids. The method integrates Kirchhoff integral formulations with thin-plate effective boundary conditions, discretizes mid-surfaces into triangular facets, and characterizes multiple reflection and transmission mechanisms through iterative pressure-field convergence. Validation against analytical solutions for a water-filled spherical shell, along with FEM simulations and experiments for a scaled BeTSSi IIB sail model, demonstrates three key outcomes:

(1) The IPA accurately reproduces interference patterns of reflected and transmitted waves, overcoming limitations of rigid and pressure-release approximations. It provides a unified framework predicting comprehensive acoustic fields, including internal cavity fields, external near-fields, and far-field scattering patterns.

(2) The IPA achieves significant computational acceleration over FEM, particularly at higher frequencies where FEM becomes prohibitively time-intensive, enabling efficient analysis of complex structures.

(3) The IPA demonstrates excellent agreement with experimental time-angle spectra, clearly resolving the primary reflection (12.6 ms) and higher-order reflection/transmission paths (>12.8 ms) in the sail model. This validates IPA's capability to characterize multiple reflection-transmission features relevant to underwater vehicle acoustics.

While the IPA framework intentionally simplifies complex elastic wave interactions and diffraction effects to maintain computational practicality for engineering design, these unmodeled phenomena contribute to localized discrepancies, particularly at low frequencies and in high-curvature regions. Future work will extend the IPA framework to elastic-dominant regimes and complex active sonar scenarios.