Study on Sound Field Properties of Parametric Array Under the Influence of Underwater Waveguide Interface Scattering Based on Non-Paraxial Model—Theory and Experiment

Abstract

This paper theoretically and experimentally studies the effect of underwater waveguide interface scattering on the nonlinear sound field characteristics of parametric array (PA) radiation. Based on the image source method, the components of the sound field in the waveguide are first analyzed. Then, a non-paraxial model is developed to account for the influence of interface scattering. This model enables accurate calculation of the wide-angle sound field. The impact of the sound source depth and the interface reflection coefficient on the distribution of the difference-frequency wave (DFW) sound field in the waveguide is studied. The interface alters the phase distribution of the DFW’s virtual source density function, thereby affecting the sound field accumulation process. Waveguide interfaces with different absorption coefficients influence the amplitude oscillation caused by interface reflection and change the sidelobe size of the DFW beam. The DFW sound field distribution is measured at three typical frequencies. Simulation and experimental results show that the attenuation of the DFW’s axial sound pressure level in the waveguide oscillates, and the DFW’s beamwidth gradually widens as the frequency decreases. The calculated results from the proposed model agree well with the measured data, with average errors along the sound axis and depth being less than 3 dB and 6 dB, respectively. This demonstrates the model’s superior applicability compared to the existing free-field model.

1. Introduction

1.1. Underwater Application of PA

PA uses the acoustic nonlinearity of the medium to generate a highly directional DFW by transmitting high-frequency large-amplitude sound waves. With the characteristics of low frequency, small size, and high directivity [1], PA has advanced applications in underwater acoustics, including ocean waveguide remote sound transmission [2,3], shallow water sub-bottom profiling [4,5], underwater acoustic material measurement [6,7], and underwater acoustic directional communication [8,9]. With the decrease in the difference frequency generated by PA, the beamwidth will gradually increase [10], making wide-angle sound field information increasingly important in PA-related applications in ocean waveguides.

1.2. Related Work

Research on the DFW sound field in waveguides began in the 1980s; Zaitsev et al. first analyzed the vertical (mode) and horizontal (angle) structure of the radiated sound field of a directional sound source [11] and then carried out relevant experimental research [12]. Li et al. investigated the effects of flow velocity and surface fluctuation on the selective excitation order of a parametric DFW sound field and proposed a Green function solution describing a DFW propagating in waveguides with low Mach number and surface fluctuation [13]. In these studies, the normal mode theory is used for analysis, and the high-frequency wave (HFW) is approximated to the quasi-direct beam to simplify the theoretical model. This simplified model is unsuitable for the diffractive PA whose HFW frequency is tens of kHz in underwater acoustics because it neglects the directional beam characteristics radiated by the sound source and thus has greater limitations. Therefore, accurately calculating the DFW sound field distribution is essential for studying the characteristics of PA-radiated sound fields in waveguides.

The theoretical modelling of the DFW sound field is more complicated than that of the linear single-frequency sound field directly radiated by the sound source. Near the sound source, the local nonlinear effects induced by the HFW are significant, typically requiring a general second-order nonlinear wave equation to model the sound propagation process [14]. In the Westervelt far-field region, the Westervelt equation, which ignores local nonlinear effects, is sufficient to describe the sound propagation process [15] accurately. However, as this is a high-order nonlinear partial differential equation with no closed-form analytical solution, numerical methods are required for its computation. When the sound pressure level (SPL) of the HFW radiated by the PA sound source is finite, the sound field governed by the Westervelt equation can be decomposed into two parts using the quasi-linear approximation and the successive approximation method. Among them, the HFW sound field can be calculated by using the Rayleigh Integral (RI) on the surface of the sound source [16], while the DFW sound field can be calculated by the virtual source volume integral on the entire observation space [10]. The resulting sound field solutions are termed quasi-linear solutions. For weak nonlinearity, the solution is accurate everywhere except near the PA source radiation surface. However, the quasi-linear solution involves a quintuple integral, resulting in significant computational costs. This makes calculations for high frequencies (tens of kHz) or large domains (hundreds of meters) highly challenging.

Subsequently, the Gaussian beam expansion (GBE) method was widely used to simplify quasi-linear solutions [17,18]. The core of the method is to approximate the vibration velocity distribution of HFW on the surface of the sound source by the sum of a set of Gaussian functions. However, the early GBE method suffers from several limitations. First, for the phased array, when the declination angle between the direction of the beam of the sound source and the sound axis is too large, the prediction accuracy on the sound axis will decline [19]. Huang et al. proposed a phase correction technique to enhance the model’s prediction accuracy [20]. However, this would cause the computational cost to increase with the increase in the array elements, and it was only effective for HFW sound fields. Second, in the wide-angle position away from the sound axis or when the frequency of the DFW is low, using the paraxial approximation to calculate the virtual source integral will lead to large errors. Cervenka, Zhuang et al. enhanced the GBE method by employing a non-paraxial approximation for the DFW virtual source integral to model the nonlinear sound field from phased PA [21,22]. However, the paraxial approximation remains for HFW calculations, leading to errors in wide-angle HFW sound field predictions. These errors significantly impact the accuracy of subsequent DFW sound field calculations.

Since the 1980s, numerous theoretical and experimental studies have explored the characteristics of PA sound sources near reflective interfaces. Muir et al. investigated the effect of pressure release surfaces on the nonlinear interaction process, demonstrating that DFW phase reversal after reflection causes interference effects [23]. Garrett studied the reflection of DFW by finite-size planar targets, but the model used was only valid in the near-field paraxial region [24]. Wang evaluated the influence of the pressure release surface on the DFW sound field of parametric source radiation from both theoretical and experimental aspects, revealing that at large grazing angles, the sound attenuation due to DFW–interface interaction approaches a constant [25]. Recently, Zhong investigated DFW reflection from PA in air [26], noting that DFW near reflective surfaces comprises both directly reflected DFW and DFW generated by reflected HFW.

1.3. Main Work of This Paper

Previous studies on DFW sound field reflection have primarily focused on a single interface, and the calculation accuracy of a wide-angle HFW sound field has not been fully addressed when calculating a nonlinear sound field. Furthermore, most studies have been conducted in an air medium, neglecting the characteristics of the water medium, particularly the waveguide environment. In this study, we develop a non-paraxial model that incorporates the scattering effects at the waveguide interface to calculate the nonlinear sound field generated by PA. First, the composition of the sound field under the influence of a waveguide is analyzed based on the image source method. Next, based on the quasi-linear solution of the second-order Westervelt equation, a non-paraxial model considering the waveguide boundary conditions is derived using only the non-paraxial approximation of the HFW sound field. The computational efficiency of the proposed model is comparable to that of the conventional GBE method in [22], and because no additional approximation is made, the sound field has better accuracy in wide-angle calculation, so that the distribution of the PA nonlinear sound field in the waveguide can be accurately obtained. The influence of the sound source’s depth and the waveguide boundary’s reflection coefficient on the distribution of the nonlinear sound field is studied. Finally, lake test results show that the proposed model has better applicability than the free-field model. The model developed in this paper can quickly obtain the PA sound field distribution in the underwater waveguide and accurately calculate the wide-angle PA sound field, which is conducive to further revealing the propagation law of nonlinear sound waves after scattering by the waveguide interface and provides theoretical support for related applications involving PA in the waveguide environment.

The paper’s structure is as follows. The second section gives the theoretical model. The third section is a simulation study of the influence factors of the nonlinear sound field. The fourth section is experimental verification. Finally, the fifth section summarizes the article.

2. Theory

In this section, the composition of the PA sound field in the waveguide is analyzed, and a non-paraxial model is established to calculate it.

2.1. Sound Field in the Waveguide

A three-dimensional rectangular coordinate system is established, as shown in Figure 1, where the $xoz$ plane is the upper interface of the waveguide, direction y is the depth direction of the waveguide, plane $y=D$ is the lower interface, and D is the total depth of the waveguide. The PA sound source (blue rectangle in Figure 1) is at depth $y=H$, its radiant plane coincides with the $xoy$ plane, and the direction of the sound axis is parallel to the z-axis. The side lengths of the PA in the x and y directions are $2a_x$ and $2a_y$, respectively. In this study, we assume that the waveguide interface is statically smooth, so the image source method can be used.

First, the HFW sound field is analyzed. In the observation area, the existence of two interfaces will make the total HFW sound field composed of the sound field radiated by the original sound source and the sound field radiated by an infinite number of image sources (grey rectangle in Figure 1), which can be expressed as

$$p_{i,sum}\left(\overrightarrow{{\mathbf{r}}_v}\right)=p_{i,initial}\left(\overrightarrow{{\mathbf{r}}_v}\right)+\sum _{\xi }{p}_{i,up}^{\xi }\left(\overrightarrow{{\mathbf{r}}_v}\right)+\sum _{\xi }{p}_{i,down}^{\xi }\left(\overrightarrow{{\mathbf{r}}_v}\right),$$

(1)

where $p_{i,initial}\left(\overrightarrow{{\mathbf{r}}_v}\right)$ is the HFW sound field radiated by the rectangular piston sound source in the free field, and $p_{i,up}^{\xi }\left(\overrightarrow{{\mathbf{r}}_v}\right)$ is the HFW sound field radiated by the $\xi$-th image source located above the $xoz$ plane. According to the geometric relationship, the coordinates of the $\xi$-th image source on the y-axis are $(0,Y_{up}^{\xi },0)$:

$$Y_{up}^{\xi }=\left\{\begin{array}{c}-\left(\right(\xi -1)D+H)\xi =1,3,5\cdots \hfill \\ -(\xi D-H)\xi =2,4,6\cdots .\hfill \end{array}\right.$$

(2)

And $p_{i,down}^{\xi }\left(\overrightarrow{{\mathbf{r}}_v}\right)$ is the HFW sound field radiated by the $\xi$-th image source below the $y=D$ plane, whose coordinates on the y-axis are $(0,Y_{down}^{\xi },0)$:

$$Y_{down}^{\xi }=\left\{\begin{array}{c}(\xi +1)D-H\xi =1,3,5\cdots \hfill \\ \xi D+H\xi =2,4,6\cdots .\hfill \end{array}\right.$$

(3)

Each additional image source is equivalent to calculating one more reflected sound ray at the surface or bottom of the water in the total sound field. The final superimposed sound fields meet the upper and lower boundary conditions of the waveguide alternately according to certain rules.

Similarly, the total DFW sound field is composed of the sound field formed by virtual sources in free space and the sound field formed by infinite image virtual sources, which can be expressed as

$$p_{d,sum}\left(\overrightarrow{\mathbf{r}}\right)=p_{d,initial}\left(\overrightarrow{\mathbf{r}}\right)+\sum _{\xi }{p}_{d,up}^{\xi }\left(\overrightarrow{\mathbf{r}}\right)+\sum _{\xi }{p}_{d,down}^{\xi }\left(\overrightarrow{\mathbf{r}}\right),$$

(4)

where $p_{d,initial}\left(\overrightarrow{\mathbf{r}}\right)$ is the DFW sound field formed by the virtual source in the free field, and $p_{d,up}^{\xi }\left(\overrightarrow{\mathbf{r}}\right)$ is the DFW sound field formed by the $\xi$-th image virtual source located above the $xoz$ plane. $p_{d,down}^{\xi }\left(\overrightarrow{\mathbf{r}}\right)$ is the DFW sound field formed by the $\xi$-th image virtual source below the $y=D$ plane. The position distribution of the image virtual source is similar to that of the HFW image source in Equations (2) and (3).

2.2. Non-Paraxial Model

The quasi-linear solution is derived from two coupled equations.

$${\nabla }^2{p}_i\left(\overrightarrow{{\mathbf{r}}_v}\right)+k_i^2{p}_i\left(\overrightarrow{{\mathbf{r}}_v}\right)=0,(i=1,2),$$

(5)

$${\nabla }^2{p}_d\left(\overrightarrow{\mathbf{r}}\right)+k_d^2{p}_d\left(\overrightarrow{\mathbf{r}}\right)=q\left(\overrightarrow{{\mathbf{r}}_v}\right),$$

(6)

where $k_i={\omega }_i/c_0+\mathrm{j}{\alpha }_i$ is the HFW complex wave number, ${\omega }_i$ is the angular frequency of the HFW, $c_0$ is the adiabatic sound velocity, $\mathrm{j}$ is the imaginary unit, ${\alpha }_i$ is the sound attenuation coefficient, $k_d$ is the DFW complex wave number, $p_i\left(\overrightarrow{{\mathbf{r}}_v}\right)$ is the HFW sound field solution, $p_d\left(\overrightarrow{\mathbf{r}}\right)$ is the DFW sound field solution, $q\left(\overrightarrow{{\mathbf{r}}_v}\right)$ is the virtual source density function, $\overrightarrow{{\mathbf{r}}_v}=(x_v,y_v,z_v)$ is the virtual source point, and $\overrightarrow{\mathbf{r}}=(x,y,z)$ is the field point.

The solution of Equation (5) can be constructed as the Rayleigh integral (RI):

$$p_i\left(\overrightarrow{{\mathbf{r}}_v}\right)=-{\displaystyle \frac{\mathrm{j}{\rho }_0{\omega }_i}{2\pi }}{\int \int }_S{v}_i(x_s,y_s){\displaystyle \frac{exp\left(\mathrm{j}{k}_i{d}_s\right)}{d_s}}\mathrm{d}{x}_s\mathrm{d}{y}_s,$$

(7)

where ${\rho }_0$ is the fluid density, $v_i\left(\overrightarrow{{\mathbf{r}}_s}\right)$ is the normal component of vibration velocity at point $\overrightarrow{{\mathbf{r}}_s}=(x_s,y_s,0)$, and $d_s$ is the distance between point $\overrightarrow{{\mathbf{r}}_s}$ on the surface of the sound source and the virtual source point $\overrightarrow{{\mathbf{r}}_v}$.

In the conventional GBE model [22], paraxial approximation is usually adopted for Equation (7). This approximation method will introduce errors at the position deviating from the z-axis, and the errors will increase with the increase in the degree of deviation. This paper considers the more general case and directly adopts the non-paraxial approximation to the HFW sound field, which can further ensure the wide-angle calculation accuracy of the PA sound field. That is, use $R_s=\sqrt{{x_v}^2+{y_v}^2+{z_v}^2}$ instead of $d_s$ in the exponential function of Equation (7). The two different approximations of $d_s$ are

$$\begin{array}{c}{d}_s=\sqrt{{\left(x_v-x_s\right)}^2+{\left(y_v-y_s\right)}^2+{\left(z_v-z_s\right)}^2}\\ \approx \left\{\begin{array}{c}\left|z_v-z_s\right|+{\displaystyle \frac{{\left(x_v-x_s\right)}^2+{\left(y_v-y_s\right)}^2}{2\left|z_v-z_s\right|}},\mathrm{conventional}\mathrm{GBE}\mathrm{model}\hfill \\ R_s+{\displaystyle \frac{{\left(x_v-x_s\right)}^2+{\left(y_v-y_s\right)}^2}{2R_s}}-{\displaystyle \frac{x_v^2+y_v^2}{2R_s}},\text{non-paraxial}\mathrm{model}\hfill \end{array}\right.\end{array}$$

(8)

The surface vibration velocity of the rectangular sound source is decomposed into:

$$\begin{array}{c}{v}_i(x_s,y_s)=V_i{\displaystyle \sum _{m=1}^N}{A}_m{e}^{\left[-B_m{x}_s^2/\phantom{B_m{x}_s^2{a_x}^2}{a_x}^2\right]}\times {\displaystyle \sum _{n=1}^N}{A}_n{e}^{\left[-B_n{y}_s^2/\phantom{B_n{y}_s^2{a_y}^2}{a_y}^2\right]},\end{array}$$

(9)

where $V_i$ is the amplitude, and $A_m$, $B_m$, and $A_n$, $B_n$ are complex coefficients, which can be obtained by numerical multi-dimensional optimization method. N is the number of Gaussian expansion coefficients. To effectively fit the distribution of vibration velocity on the surface of a rectangular sound source, generally, $N=10$ can meet the requirements [19].

The following dimensionless variables are introduced to simplify the calculation:

$$\left\{\begin{array}{c}\widetilde{x_s}=x_s/\phantom{x_s{a}_x}{a}_x,\widetilde{y_s}=y_s/\phantom{y_s{a}_y}{a}_y,\widetilde{x_v}=x_v/\phantom{x_v{a}_x}{a}_x,\widetilde{y_v}=y_v/\phantom{y_v{a}_y}{a}_y\hfill \\ \widetilde{z_{v,\eta ,i}}=z_v/R_{\eta ,i},\widetilde{R_{s,\eta ,i}}=R_s/R_{\eta ,i},R_{\eta ,i}=k_i{a_{\eta }}^2/2\hfill \\ (\eta =x,y;i=1,2).\hfill \end{array}\right.$$

(10)

By substituting Equations (8)–(10) into Equation (7), the HFW sound field solution under non-paraxial approximation can be obtained through derivation:

$$\begin{array}{c}{p}_i\left(\overrightarrow{{\mathbf{r}}_v}\right)=P_i{e}^{\left\{\mathrm{j}{k}_i\left[R_s-\left({x_v}^2+{y_v}^2\right)/\phantom{\left({x_v}^2+{y_v}^2\right)\left(2R_s\right)}\left(2R_s\right)\right]\right\}}\times \\ {\displaystyle \sum _{m_i=1}^N}\left[A_{m_i}\times G\left(\widetilde{x_v},\widetilde{R_{s,x,i}},B_{m_i}\right)\right]\times \\ {\displaystyle \sum _{n_i=1}^N}\left[A_{n_i}\times G\left(\widetilde{y_v},\widetilde{R_{s,y,i}},B_{n_i}\right)\right].\end{array}$$

(11)

When considering the influence of waveguide, it is necessary to consider the boundary reflection condition of the sound field radiated by each image source in the image process. Taking the first image source under the $y=D$ plane as an example, the HFW sound field formed by it is

$$\begin{array}{c}{p}_{i,down}^1\left(\overrightarrow{{\mathbf{r}}_v}\right)=-{\displaystyle \frac{\mathrm{j}{\rho }_0{\omega }_i}{2\pi }}{\int \int }_S[R_{i,down}\left(w_i\right)v_i(x_s,y_s){\displaystyle \frac{e^{\left(\mathrm{j}{k}_i{d}_{i,down}^1\right)}}{d_{i,down}^1}}]\mathrm{d}{x}_s\mathrm{d}{y}_s,\end{array}$$

(12)

where $R_{i,down}\left(w_i\right)$ is the reflection coefficient of the lower interface, and the distance $d_{i,down}^1$ is

$$d_{i,down}^1=\sqrt{{(x_v-x_s)}^2+{(y_v-y_s-Y_{down}^1)}^2+{(z_v-z_s)}^2}.$$

(13)

Equations (12) and each of the terms of (1) can be derived in a form similar to Equation (11), from which a non-paraxial model of the HFW sound field generated by PA in the waveguide is obtained.

The DFW sound pressure radiated by PA can be thought of as a superposition of the sound pressure produced by an infinite virtual source in the medium

$$p_d\left(\overrightarrow{\mathbf{r}}\right)=-{\displaystyle \frac{\mathrm{j}{\rho }_0{\omega }_d}{4\pi }}{\int \int \int }_{V}q\left(\overrightarrow{{\mathbf{r}}_v}\right){\displaystyle \frac{exp\left(\mathrm{j}{k}_d{d}_v\right)}{d_v}}\mathrm{d}{x}_v\mathrm{d}{y}_v\mathrm{d}{\mathrm{z}}_v,$$

(14)

where $d_v$ is the distance between the virtual source point $\overrightarrow{{\mathbf{r}}_v}$ and the spatial field point $\overrightarrow{\mathbf{r}}$. The virtual source density function is:

$$q\left(\overrightarrow{{\mathbf{r}}_v}\right)=\beta {\omega }_d^2{p}_1\left(\overrightarrow{{\mathbf{r}}_v}\right)p_2^{*}\left(\overrightarrow{{\mathbf{r}}_v}\right)/{\rho }_0{c}_0^4.$$

(15)

By substituting Equations (11) and (15) into Equation (14), the DFW sound field solution under non-paraxial approximation can be obtained through derivation:

$$\begin{array}{c}{p}_d\left(\overrightarrow{\mathbf{r}}\right)=-{\displaystyle \frac{\beta {\omega }_d^2{P}_1{P}_2}{4\pi {\rho }_0{c}_0^4}}\underset{V}{\int \int \int }{e}^{\left\{\mathrm{j}{k}_d\left[R_s-\left({x_v}^2+{y_v}^2\right)/\phantom{\left({x_v}^2+{y_v}^2\right)\left(2R_s\right)}\left(2R_s\right)\right]\right\}}{\displaystyle \frac{{\mathrm{e}}^{\mathrm{j}{k}_d{d}_v}}{d_v}}\times \hfill \\ {\displaystyle \sum _{{\displaystyle \genfrac{}{}{0pt}{}{m_1=1\hfill }{m_2=1\hfill }}}^N}{A}_{m_1}G(\widetilde{x_v},\widetilde{R_{s,x,1}},B_{m_1}){\left[A_{m_2}G(\widetilde{x_v},\widetilde{R_{s,x,2}},B_{m_2})\right]}^{*}\times \hfill \\ {\displaystyle \sum _{{\displaystyle \genfrac{}{}{0pt}{}{n_1=1\hfill }{n_2=1\hfill }}}^N}{A}_{n_1}G(\widetilde{y_v},\widetilde{R_{s,y,1}},B_{n_1}){\left[A_{n_2}G(\widetilde{y_v},\widetilde{R_{s,y,2}},B_{n_2})\right]}^{*}\mathrm{d}\overrightarrow{{\mathbf{r}}_v},\hfill \end{array}$$

(16)

where $G(\zeta ,\sigma ,\gamma )=exp\left[-\gamma {\zeta }^2/\phantom{\gamma {\zeta }^2\left(1+\mathrm{j}\gamma \sigma \right)}\left(1+\mathrm{j}\gamma \sigma \right)\right]/\sqrt{1+\mathrm{j}\gamma \sigma }$.

When considering the influence of waveguide, taking the first image virtual source under the $y=D$ plane as an example, the DFW sound field formed by it is

$$\begin{array}{c}{p}_{d,down}^1\left(\overrightarrow{\mathbf{r}}\right)=-{\displaystyle \frac{\mathrm{j}{\rho }_0{\omega }_d}{4\pi }}\underset{V}{\int \int \int }[R_{d,down}\left(w_d\right)q_d\left(\overrightarrow{{\mathbf{r}}_v}\right){\displaystyle \frac{e^{\left(\mathrm{j}{k}_d{d}_{v,down}^1\right)}}{d_{v,down}^1}}]\mathrm{d}V,\end{array}$$

(17)

where $R_{d,down}\left(w_d\right)$ is the reflection coefficient of the lower interface, and the distance $d_{v,down}^1$ is

$$\begin{array}{c}{d}_{v,down}^1=\sqrt{{(x-x_v)}^2+{(y-y_v-Y_{down}^1)}^2+{(z-z_v)}^2}.\end{array}$$

(18)

Equations (17) and each of the terms of (4) can be derived in a form similar to Equation (16), from which a non-paraxial model of the DFW sound field generated by PA in the waveguide is obtained.

By comparing Equations (7) and (11), as well as Equations (14) and (16), it can be seen that the expression of the HFW sound field is transformed from a double integral into a summation form, and the solution of the DFW sound field is also simplified from the original five-dimensional integral to a triple integral, which will provide convenience for the subsequent solution of the DFW sound field. The computational cost of this model is equivalent to that of the traditional GBE method, and the wide-angle PA sound field can be accurately calculated by breaking the limitation of paraxial approximation.

3. Simulations and Discussion

In this section, the non-paraxial model proposed in Section 2.2 is verified, and the sound field in the waveguide is simulated by using this model.

3.1. Verification of Non-Paraxial Model

The non-paraxial model proposed in Section 2.2 adopts strict virtual source integration in calculating the DFW sound field without any additional approximation. The calculation accuracy of the HFW sound field will directly reflect the calculation accuracy of the final DFW sound field. Therefore, the HFW sound field is studied here. The result obtained by the Rayleigh integral model (RI, Equation (7)) is taken as the exact solution, and the conventional GBE model [22] is used as a comparison to illustrate the applicability of the proposed model.

The simulation parameters are set as follows. In the three-dimensional free-field homogeneous water medium, the acoustic parameters of the medium do not vary with the spatial position. The sound velocity of the water medium is $c_0=1482\mathrm{m}/\mathrm{s}$, the density is ${\rho }_0=998\mathrm{kg}/{\mathrm{m}}^3$, the nonlinear parameter is $\beta =3.5$, the reference sound pressure in the water medium is 1 μPa, and the attenuation is calculated according to the empirical formula of water absorption [27]. The vibration velocity of the PA sound source surface is uniformly distributed and the peak sound pressure on the surface of the sound source is 100 $\mathrm{kPa}$. The side length of the sound source is $L_x=0.072\mathrm{m}$, $L_y=0.1\mathrm{m}$, and $L_c=(L_x+L_y)/2$ is defined. The HFW frequencies are $f_1=40.1\mathrm{kHz}$ and $f_2=39.9\mathrm{kHz}$, respectively, and the DFW frequency is $f_d=200\mathrm{Hz}$. The HFW wavelengths are ${\lambda }_1=0.0370\mathrm{m}$ and ${\lambda }_2=0.0371\mathrm{m}$, respectively. The Cartesian rectangular coordinate system $O-xyz$ is established with the geometric center of the PA radiation surface as the origin, the sound axis coincides with the z axis, and the spatial range of the sound field is $x\in [-10\mathrm{m},10\mathrm{m}]$, $y\in [-10\mathrm{m},10\mathrm{m}]$, $z\in [0\mathrm{m},20\mathrm{m}]$. The geometric center of the PA radiation surface is located at $(0,0,0)$. The computing platform is based on AMD Ryzen Threadripper PRO 5955WX, which has 16 cores, 4.00 GHz, and 128 GB RAM.

In the simulation, the HFW with frequency $f_1$ is taken as the research object. The numerical solution of the RI model is obtained by gradually reducing the integral step of the sound source surface. The result converges when the sound pressure level (SPL) error is less than 0.1 dB. Finally, the integral step of the RI method is ${\lambda }_1/6$. The sound pressure level (SPL) distribution calculated by the three models is shown in Figure 2, where the $yoz$ cross-section is taken for analysis. In the free field, the HFW directly radiated by a finite-size piston sound source will form a directional beam with a sidelobe in space, consistent with the results in Figure 2. By comparing Figure 2a,b, it can be found that the HFW sound field obtained by the conventional GBE method has apparent errors in the wide-angle position. By comparing Figure 2a,c, it can be found that the sound field calculated by the non-paraxial model is closer to the exact solution than that calculated by the conventional GBE model, especially at the wide-angle position away from the sound axis.

We calculate the axial SPL, sound vertical directivities at 10 m from the source, and the errors between them and the exact solution, as shown in Figure 3, where the SPL error is defined as ${\Delta }_{\mathrm{SPL}}=\left|20lg\left(\left|p^{\prime }\right|/\phantom{\left|p^{\prime }\right|\left|p\right|}\left|p\right|\right)\right|$. Figure 3a shows that the HFW follows the spherical wave attenuation rule, and Figure 3b shows that it has sharp directivity but also has sidelobe. According to the analysis in Section 2, the calculation results of the non-paraxial model and the conventional GBE model on the sound axis are the same. In Figure 3b, the directivity calculated by the non-paraxial model is very consistent with the exact RI solution within the range of $\pm {20}^{\circ }$ and is basically consistent with the larger angle, while the results of the conventional GBE model have errors within $\pm {20}^{\circ }$, and the errors gradually become larger as the angle increases.

In Figure 3c, the “effective region” (area restricted by white dotted lines in Figure 3c,d) of the conventional GBE model has only a range of $\pm {15}^{\circ }$, while in Figure 3d, the “effective region” of the non-paraxial model has a range of $\pm {20}^{\circ }$. It is obviously observed that the accuracy of the non-paraxial model is better than that of the conventional GBE model at a large angle. It can be seen that the conventional GBE model is only suitable for the calculation of the HFW sound field in the paraxial region (usually refers to $\pm {15}^{\circ }$), while the non-paraxial model is suitable for the calculation of the wide-angle HFW sound field. Based on the analysis in Section 2.2, the discrepancy in approximation accuracy arises from the different reference points chosen by the two approximation methods during the second-order Taylor expansion of $d_s$.

Table 1. Calculation time of different models.

	${\mathit{f}}_{\mathit{c}}$ = 40 kHz ${\mathit{L}}_{\mathit{c}}$ = 0.086 m $\mathit{ka}$ = 14.6	${\mathit{f}}_{\mathit{c}}$ = 60 kHz ${\mathit{L}}_{\mathit{c}}$ = 0.086 m $\mathit{ka}$ = 21.9	${\mathit{f}}_{\mathit{c}}$ = 40 kHz ${\mathit{L}}_{\mathit{c}}$ = 0.172 m $\mathit{ka}$ = 29.2
RI	75 s	591 s	317 s
non-paraxial	28 s	97 s	29 s
conventional GBE	27 s	95 s	28 s

shows the calculation time required by the three models under different $ka$ values. The non-paraxial model can effectively reduce the calculation time cost compared with the RI model. Therefore, the subsequent calculation in this paper will be based on the non-paraxial model.

3.2. Study on PA Sound Field Characteristics in Waveguide

3.2.1. The Influence of Waveguide on the Vertical Directivity of Sound Source

According to the proposed model, the distribution of the DFW field generated by the PA in the waveguide with a frequency of 200 Hz is calculated. Where the sound source and environmental parameters are consistent with those in Section 3.1, the coordinates of the PA sound source are $(0\mathrm{m},50\mathrm{m},0\mathrm{m})$. The spatial range of sound field observation is as follows: $x\in [-50\mathrm{m},50\mathrm{m}]$, $y\in [0\mathrm{m},100\mathrm{m}]$, $z\in [0\mathrm{m},100\mathrm{m}]$. For ideal underwater waveguides, the upper interface of the waveguide is the Dirichlet boundary, and the sound reflection coefficient is $-1$. The lower interface is the Neumann boundary, and the sound reflection coefficient is 1.

The linear low-frequency sound field directly generated by the sound source with a frequency of 200 Hz is compared, and the sound vertical directivity at the distance of 50 m from the sound source in the free field and waveguide is calculated. The results are shown in Figure 4. For the linear sound field directly generated by the sound source, Figure 4c shows that the low-frequency sound wave is almost fully directional and uniform in the free field but will show a strong interference phenomenon in the waveguide, which will make its vertical directivity uneven. Figure 4a,b show the amplitude and phase distribution of the linear low-frequency sound field in the waveguide, which appears to be disorganized due to multiple reflections at the interface. For the DFW sound field generated by the sound source, the DFW has a sharp vertical directivity and no sidelobe in the free field, and the energy of the sound field is mainly concentrated near the sound axis. The waveguide will also cause interference in the sound field, and Figure 4f shows that the influence brought by such interference will introduce sidelobe in the wide-angle position. By comparing Figure 4a,b,d,e, it can be found that the influence of the waveguide on the DFW sound field is less significant than that of the linear sound field.

3.2.2. The Influence of Sound Source Depth on Sound Field Distribution

Because of the sharp directivity of the DFW beam, when the distance between the sound source and the waveguide boundary is far, the beam does not fully interact with the boundary, resulting in the sound field near the sound axis not changing significantly and the vertical directivity of the sound field is not much different from that of the free field. In order to further study the influence of the waveguide boundary, the distribution of the DFW sound fields generated by PA sound sources close to the upper and lower interfaces of the waveguide is calculated, respectively, and the results are shown in Figure 5. By comparing Figure 5a,c,d,f, it can be found that when the sound source is close to the waveguide boundary, the DFW beam will be significantly affected by the boundary, and whether it is an “absolutely soft” boundary or an “absolutely hard” boundary, the beam will show sidelobe in the region near the boundary, and cause the vertical directivity to lose symmetry. Figure 5b,e show that the spatial phase of the DFW radiated by the PA sound source at different depths is more obviously affected by the interface away from the sound source side and thus produces chaotic oscillating fringes.

To study the mechanism of sidelobe generation caused by the waveguide boundary, the amplitude and phase of the integrand function of the virtual source integral in Equation (14) are simulated when the sound source is close to the upper and lower interfaces of the waveguide, respectively, and the free field is compared. The results are shown in Figure 6. When the sound source is located at 10 m and 90 m, the coordinates of any field point in the observation area are selected as $(0\mathrm{m},10\mathrm{m},50\mathrm{m})$ and $(0\mathrm{m},90\mathrm{m},50\mathrm{m})$, respectively.

As shown in Figure 6a,c,e,g, the amplitude distribution of the integrand function has a greater weight in the center of the sound source radiation plane and the field point, and Figure 6b,f show that the phase of the integrand function is relatively uniform in the area near the “source-field point line” connected between the field point and the center of the sound source radiation plane. At the same time, violent oscillations will occur outside the region, indicating that the “source-field point line” region is the key region for virtual sound source accumulation, and any small change in this region will significantly impact the final DFW sound field distribution. Figure 6d,h show that waveguide boundaries can disrupt the phase uniformity of the region, especially near the sound source. The phenomenon that the waveguide boundary causes the sidelobe of the DFW beam can be attributed to the fact that the influence of the waveguide boundary changes the phase distribution of the integrated function, so the existence of the waveguide boundary will affect the accumulation of virtual sources in the observation region, and thus affect the final distribution of the nonlinear sound field.

3.2.3. The Influence of Boundary Reflection Coefficient on Sound Field Distribution

The upper and lower interfaces of the actual underwater waveguide are not “absolutely soft” and “absolutely hard”. The following will consider the impedance characteristics of the actual medium, recalculate the sound reflection coefficient of the interface, and substitute it into the calculation of the PA sound field in the waveguide. The reflection coefficient of the plane sound wave when it is vertically incident on the interface of any two layers of media is

$$r_p={\displaystyle \frac{Z_2-Z_1}{Z_2+Z_1}},$$

(19)

where $Z_1$ and $Z_2$ are the characteristic impedances of the two dielectric layers, respectively.

Because the characteristic impedance difference between the water medium and air medium is too large, the variation range of the sound reflection coefficient of the water surface is very small in the actual situation. This paper focuses on the effect of the lower boundary of the waveguide with different materials on the DFW sound field.

Assuming that the sound velocity of air medium $c_a=343\mathrm{m}/\mathrm{s}$, the density ${\rho }_a=1.21\mathrm{kg}/{\mathrm{m}}^3$, and the parameters of the water medium are the same as in Section 3.1, the acoustic parameters of the riverbed medium (made of basalt or silt) are derived from [28], and the depth of the sound source is 90 m. The reflection coefficient of the air–water boundary is $-0.9995$, and that of the water–riverbed boundary is 0.8086 (basalt) and 0.2819 (silt), respectively. The axial SPL distribution at the impedance interface and the vertical directivity of the DFW sound field 10 m away from the sound source are obtained by simulation. At the same time, the results of the free-field case and ideal absolutely hard boundary are also shown in Figure 7 as a reference.

By comparing the sound field results under different interface conditions, it can be found that with the increase in the characteristic impedance of the medium, the sound absorption coefficient of the lower interface gradually decreases, and it will approach the characteristic of the “absolutely hard” boundary. As a result, the stronger the oscillation effect caused by the boundary reflection of the sound field, the larger the sidelobe of the DFW beam.

4. Experiment and Analysis

4.1. Experimental Setup

We conducted experiments in the waters of Qiandao Lake, as shown in Figure 8, to further verify the validity of the proposed model. The PA sound source used in the experiment is a custom rectangular piston transducer. The receiving system consists of a nine-element vertical linear array of RESON TC4013 hydrophones and a PULSE analyzer (Brüel & Kjær Type 3052, Naerum, Denmark) placed on a movable boat, where the element of the vertical array is 0.75 m apart. The instrument connection is shown in Figure 9.

Figure 8 marks the location of the measuring points selected in the experiment. The distance between the 13 measuring points and the sound source from near to far in the direction of the sound axis is shown in

Table 2. Distance between each measuring point and sound source.

number	1	2	3	4	5	6	7
distance $\left(\mathrm{m}\right)$	5	55	116	163	223	291	343
number	8	9	10	11	12	13
distance $\left(\mathrm{m}\right)$	383	440	517	624	877	1000

. In the experiment, the sound source is placed at the water depth of 5 m below platform A to transmit a dual-frequency modulated signal to generate a nonlinear sound field in space, and then the vertical array is placed at each measuring point to the water depth of 5 m to receive the sound signal (the depth of the central array element is 5 m).

In addition, the RBRduet³ T.D logger measures the water depth and temperature at each measuring point the boat moves to. Figure 10 shows the data recorded by the logger. The average water depth of the first nine measuring points is about 22 m, and the average water temperature is about 20 °C. From the 10th measuring point, the depth and water temperature fluctuate considerably.

In the experiment, the sound source transmits the dual-frequency continuous wave signal at the center frequency of 88 kHz, and three kinds of difference frequencies are tested in total. The specific emission parameters are shown in

Table 3. Emission parameters of sound source.

${\mathit{f}}_1\left(\mathbf{kHz}\right)$	${\mathit{f}}_2\left(\mathbf{kHz}\right)$	${\mathit{f}}_{\mathit{d}}\left(\mathbf{kHz}\right)$
88	95.98	7.98
88	91.98	3.98
88	89.98	1.98

.

4.2. Experimental Results

The sound field pressure data are obtained by analyzing and processing the collected signal. The DFW sound field in water with a width of 60 m, a depth of 22 m and a length of 1000 m is simulated according to the proposed model combined with experimental emission parameters and environmental parameters. Firstly, the axial SPL distribution of the DFW sound field between simulation and measurement is compared. The error between the simulation results of the proposed model and the experimental data is calculated, and the error shadow distribution is drawn (Figure 11d,e,f represent the experimental data by blue dots, and the red area is the error shadow between the calculated results of the proposed model and the experimental data). At the same time, the results calculated by the free-field model are taken as a reference, and the results are shown in Figure 11.

As shown in Figure 11a–c, for three different frequencies, the axial SPL of DFW obtained from the free-field model exhibits a smooth attenuation trend with distance. In contrast, the axial SPL of DFW in the waveguide, calculated using the proposed model and measured data, demonstrates an oscillatory attenuation pattern. Based on the findings in Section 3.2, this phenomenon can be attributed to the interference fringes generated by the interaction between the beam and the waveguide interface. At the first nine measuring points, the results obtained by the proposed model are significantly closer to the measured data than those obtained by the free-field model. In addition, by observing the error shadow distribution between the results obtained by the proposed model and the measured data in Figure 11d–f, it can be found that the error is relatively small for the three typical frequencies at the first nine measuring points. But from the tenth point on, the error starts to increase. According to Figure 10, the analysis may be caused by sudden changes in water depth and water temperature at the 10th position and drastic changes in water depth and water temperature at subsequent measuring points. In addition, it can be seen from the location distribution of measuring points marked in Figure 8 that the islands block the 12th and 13th measuring points. Moreover, it is increasingly challenging to align the sound axis of the receiving array with that of the sound source (due to measuring points 12 and 13 gradually deviating from the sound axis). These factors lead to a significant increase in the error of the results at a certain distance. Based on the data from the first nine effective measuring points, the average error of the proposed model and the free-field model relative to the measured data is shown in

Table 4. Average error (dB) of SPL for the first 9 measuring points on sound axis.

	${\mathit{f}}_{\mathit{d}}=7.98\mathbf{kHz}$	${\mathit{f}}_{\mathit{d}}=3.98\mathbf{kHz}$	${\mathit{f}}_{\mathit{d}}=1.98\mathbf{kHz}$
proposed model	1.9	2.5	2.7
free-field model	4.8	2.6	20.1

. For three different frequencies, the average error of the proposed model and the measured data is less than 3 dB, which has better applicability than the free-field model, and the free-field model cannot reflect the oscillation characteristics of the sound field changes.

Accurately testing the directivity of a sound source in a non-laboratory environment is challenging, so we obtain the sound field distribution in the depth direction of the DFW by using the signals received by different elements of a vertical array to varying depths of the same measuring point and use it to replace the vertical directivity. The sound field distribution of DFW at the measuring point of 5 m was selected to be studied. At this time, the test angle covering the sound field was about ±31°. The average error between the proposed model and the measured data was calculated, and the results in the free field were taken as a reference, as shown in Figure 12 and

Table 5. Average error (dB) of SPL in depth direction at 5m measuring point.

	${\mathit{f}}_{\mathit{d}}=7.98\mathbf{kHz}$	${\mathit{f}}_{\mathit{d}}=3.98\mathbf{kHz}$	${\mathit{f}}_{\mathit{d}}=1.98\mathbf{kHz}$
proposed model	5.1	4.5	2.4
free-field model	5.8	4.8	2.9

. In Figure 12a–c, it is evident that the beamwidth of the DFW gradually increases as the frequency decreases. As analyzed in Section 3.2, when the PA sound source is positioned close to the waveguide interface, scattering at the interface significantly affects the accumulation of DFW virtual sources in the region. Additionally, during experimental testing, fluctuations in the water surface cause the vertical array to sway, introducing further measurement uncertainties. Figure 12 and Table 5 show that with the increase in the frequency of the DFW, the error between the calculated results obtained by using the proposed method and the measured data gradually increases, which may be because the beamwidth of the DFW with higher frequency is sharper, which makes it more difficult for the receiving array to align the sound axis. Nevertheless, the calculated results obtained by using the proposed model are basically consistent with the measured data, and the average error between the proposed model and the measured data is less than 6 dB.

5. Conclusions and Discussion

This paper presents a theoretical and experimental study of the DFW sound field generated by PA in underwater waveguides. A non-paraxial model that accounts for the effects of waveguide boundary scattering is developed. The main findings are summarized as follows:

• The two interfaces of the waveguide result in the HFW (DFW) sound field being a combination of the sound field generated by the actual source (or virtual source) in free space and the sound fields generated by infinite image sources (or image virtual sources). When calculating the sound field radiated by each image source (or image virtual source), the boundary reflection conditions imposed during the imaging process must be taken into account.
• The characteristics of linear low-frequency sound field and nonlinear DFW sound field generated by the sound source in the waveguide are investigated, and the effects of the depth of the sound source and the boundary reflection coefficient of the waveguide on the distribution of the nonlinear sound field are analyzed. The waveguide boundary influences the virtual source accumulation process by disrupting the phase uniformity of the virtual source integral, thereby altering the final distribution of the DFW sound field.
• The sound absorption coefficients of waveguide interfaces vary with their characteristic impedances. These variations influence the intensity of the amplitude oscillations caused by interface reflections and alter the sidelobe size of the DFW beam’s directivity.
• The spatial distribution of DFW sound fields at different frequencies is measured experimentally. In the direction of the sound axis, the average error between the proposed model and the measured data is less than 3 dB; in the direction of depth, the average error between the proposed model and the measured data is less than 6 dB, effectively verifying the applicability of the proposed model.

The model developed in this study enables rapid determination of the DFW sound field distribution generated by PA in underwater waveguide environments and provides accurate calculations of wide-angle PA sound fields. This facilitates a deeper understanding of the propagation laws of nonlinear sound waves following scattering at the waveguide interface. It is important to note that the proposed model relies on the validity of the quasi-linear method. However, the conversion efficiency of PA has always attracted much attention. In the practical application of PA, it is always hoped to improve the DFW amplitude as much as possible. Although the DFW amplitude will increase with the SPL of HFW within a specific range, the nonlinearity will become strong when the SPL of HFW is too large, and the waveform will be seriously distorted. Additionally, shock waves may form during propagation [29]. At this time, the quasi-linear method is no longer applicable, and the nonlinear sound field distribution must be solved by other methods [30,31]. The Goldberg number [32] is often used to determine whether nonlinear distortion has occurred. In future studies, we will investigate the distortion of nonlinear sound waves generated by PA in waveguides. Another noteworthy problem is that the image source method employed in this paper is usually used to solve the problem of static smooth interface reflection and propagation. When the interface undergoes depth-wise motion or displacement, the image source experiences a frequency shift due to the Doppler effect, thereby altering the observed frequency. In subsequent work, we will consider using the modified image source method [33] to study this problem. If a rough interface is involved, the image source method is no longer applicable because there is no one-to-one correspondence between the incident and reflected sound rays when reflections along the rough interface are considered. In future studies, we will explore the use of the finite-difference method [34] to analyze PA sound field distributions in dynamic, rough interfaces and inhomogeneous waveguide media.