A Chebyshev–Tau Spectral Method for Coupled Modes of Underwater Sound Propagation in Range-Dependent Ocean Environments

Houwang Tu [email protected] College of Computer, National University of Defense Technology, Changsha, 410073, China Yongxian Wang [email protected] College of Meteorology and Oceanography, National University of Defense Technology, Changsha, 410073, China Chunmei Yang [email protected] First Institute of Oceanography and Key Laboratory of Marine Science and Numerical Modeling, Ministry of Natural Resources, Qingdao, 266061, China Wei Liu College of Meteorology and Oceanography, National University of Defense Technology, Changsha, 410073, China Wenbin Xiao College of Meteorology and Oceanography, National University of Defense Technology, Changsha, 410073, China Xiaodong Wang College of Computer, National University of Defense Technology, Changsha, 410073, China

Abstract

ABSTRACT The stepwise coupled-mode model is a classic approach for solving range-dependent sound propagation problems. Existing coupled-mode programs have disadvantages such as high computational cost, weak adaptability to complex ocean environments and numerical instability. In this paper, a new algorithm is designed that uses an improved range normalization and global matrix approach to address range dependence in ocean environments. Due to its high accuracy in solving differential equations, the spectral method has recently been applied to range-independent normal modes and has achieved remarkable results. This algorithm uses the Chebyshev–Tau spectral method to solve for the eigenmodes in the range-independent segments. The main steps of the algorithm are parallelized, so OpenMP multithreading technology is also applied for further acceleration. Based on this algorithm, an efficient program is developed, and numerical simulations verify that this algorithm is reliable, accurate and capable. Compared with the existing coupled-mode programs, the newly developed program is more stable and efficient at comparable accuracies and can solve waveguides in more complex and realistic ocean environments.

Keywords: Spectral method; coupled modes; range dependent; underwater acoustics; computational ocean acoustics.

I Introduction

The numerical sound field of a range-dependent waveguide is a research hot spot in computational ocean acoustics. At present, techniques for solving range-dependent acoustic propagation problems include coupled modes, adiabatic modes, rays Jensen2011; Etter2018, the parabolic approximation RAM and direct solutions to the Helmholtz equation using finite difference Liuw2021 or finite element methods Murphy1988a; Murphy1988b; Murphy1989; Murphy1996. Each method or model has its own advantages and disadvantages. Coupled-mode theory is a classic model to solve sound propagation problems in range-dependent ocean environments. It is often used to provide benchmark solutions to test the reliability of other numerical models because of its high accuracy.

The classic normal mode theory proposed by Pekeris Pekeris1948 provides solutions suitable only for range-independent acoustic waveguides and is powerless for range-dependent problems. The theory of coupled modes was proposed by Pierce Pierce1954 and Miller Miller1954 in 1954; they asserted that energy is exchanged between normal modes in a horizontally changing waveguide. Subsequently, Rutherford and Hawker Rutherford1981 noted that Pierce and Miller’s use of vertical derivative operators to replace normal derivative operators resulted in nonconservation of energy in sloping terrains; consequently, they proposed a first-order modification to coupled-mode theory to maintain first-order conservation of energy on slopes; Fawcett provided a full, analytically exact evaluation of these same terms Fawcett1992. In 1983, Evans Evans1983 proposed the idea of using a stair-step geometry to discretize sloping terra, where each step was considered a flat segment. In combination with boundary conditions, the propagator matrix between the coupling coefficients of each segment can be obtained, and the coupling coefficients of the segments can be obtained by considering radiation conditions. The acoustic field solution of each segment contains both the forward scattering mode, which exponentially decays with increasing range, and the backward scattering mode, which exponentially grows with increasing range. When considering leaky modes, the traditional superposition method suffers from numerical instability. In 1985, Mattheij Mattheij1985 proposed a decoupling matrix algorithm to solve the two-point boundary value problem. Soon after, Evans Evans1986 applied this decoupling algorithm to stepwise coupled modes, successfully resolved the numerical instability caused by leaky modes, and developed the numerical program COUPLE. The latest version, COUPLE07 Couple, can accurately calculate the fully elliptic two-way solution of the Helmholtz equation, which is considered to be an outstanding representative of coupled modes and has been widely used for many years to provide accurate solutions for numerical experiments.

However, Luo et al. Luowy2012a; Luowy2012b; Luowy2012c; Luowy2012d and Yang et al. Yangcm2012; Yangcm2015a reported that COUPLE exhibited numerical instability due to unreasonable normalized range solutions. In solving for the range-independent normal modes, COUPLE employs the Galerkin method, which forms a generalized eigenvalue problem of symmetric matrices $\mathbf{A}$ and $\mathbf{B}$ Couple in each segment:


	$\displaystyle\mathbf{Au}=\lambda\mathbf{Bu}$		(1a)
	$\displaystyle A_{n,m}=\int_{0}^{H}\frac{k^{2}(z)\phi_{n}(z)\phi_{m}(z)}{\rho(z)}\mathrm{d}z-\int_{0}^{H}\frac{1}{\rho(z)}\frac{\mathrm{d}\phi_{n}(z)}{\mathrm{d}z}\frac{\mathrm{d}\phi_{m}(z)}{\mathrm{d}z}\mathrm{d}z$		(1b)
	$\displaystyle B_{n,m}=\int_{0}^{H}\frac{\phi_{n}(z)\phi_{m}(z)}{\rho(z)}\mathrm{d}z$		(1c)

where $\phi(x)$ are the basis/weight functions in the Galerkin method. Although the matrices $\mathbf{A}$ and $\mathbf{B}$ are both formally symmetrical (symmetry means that such a generalized eigenvalue problem is efficient to solve), the elements in matrices $\mathbf{A}$ and $\mathbf{B}$ must be individually obtained through numerical quadrature, which requires many calculations. In addition, the Galerkin method must construct basis functions that satisfy the boundary conditions in each segment, which imposes considerable computational cost. Furthermore, COUPLE considers only two layers of media, which is a limitation in many cases. For example, for the lower boundary of the acoustic half-space, the bottom sediment of COUPLE needs to be set as an absorbing layer, which precludes flexibility for complicated waveguides. The KRAKEN program based on the finite difference method has good flexibility in solving for range-independent normal modes, but it can calculate only one-way coupled modes, and the stability of the coupling is often unsatisfactory Kraken2001.

In recent years, many studies have begun to address acoustic propagation problems by applying more accurate spectral methods Tuhw2020a; Tuhw2021a; Tuhw2021b; Wangyx2021a; SMPE; Tuhw2021c; Wangyx2021b. In 1993, Dzieciuch Dzieciuch1993; aw first used the Chebyshev–Tau spectral method to solve for the normal modes of the water column. Evans rimLG in 2016 devised a Legendre–Galerkin spectral method to solve the problem of acoustic propagation in a two-layer ocean environment that contained bottom sediment. In 2020, Tu et al. Tuhw2020a; Tuhw2021a used the Chebyshev–Tau spectral method to more efficiently solve this problem. Numerical experiments have shown that the NM-CT program based on the Chebyshev–Tau spectral method NM-CT is faster than the rimLG program rimLG based on the Legendre–Galerkin spectral method and more accurate than the classic finite difference method Kraken2001. Recently, Sabatini et al. Sabatini2019 and Tu et al. Tuhw2021c; MultiLC used the Chebyshev collocation method and Legendre collocation method, respectively, to solve the problem of acoustic propagation in multilayer media. Existing studies have shown that spectral methods can solve underwater acoustic propagation problems with high accuracy. However, the current programs aw; rimLG; NM-CT; MultiLC based on spectral methods can provide solutions for only range-independent acoustic waveguides. The present article combines stepwise coupled modes with the Chebyshev–Tau spectral method to develop a new algorithm that can efficiently provide solutions for range-dependent acoustical waveguides. Compared with the existing program-based coupled modes, the capability and computational efficiency of the algorithm proposed in this paper are greatly improved while maintaining the same accuracy.

II Physical Model

II.1 Range-independent normal modes

We consider a two-dimensional point source acoustic field in a cylindrical axisymmetric environment, where the angular frequency of the acoustic source is $\omega$ and the simple harmonic point source is located at $r=0$ with $z=z_{\mathrm{s}}$ . Let the acoustic pressure be $p=p(r,z)$ , and omit the time factor $\exp(-\mathrm{i}\omega t)$ . The acoustic governing equation (Helmholtz equation) can be written as Jensen2011:

\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial p}{\partial r}\right)+\rho(z)\frac{\partial}{\partial z}\left(\frac{1}{\rho(z)}\frac{\partial p}{\partial z}\right)+\frac{\omega^{2}}{c^{2}(z)}p=-\frac{\delta(r)\delta(z-z_{\mathrm{s}})}{2\pi r}

(2)

where $\omega=2\pi f$ , $f$ is the frequency of the sound source, and $c(z)$ and $\rho(z)$ are the sound speed and density profiles, respectively.

Using the technique of the separation of variables Pekeris1948, the acoustic pressure can be decomposed into:

p(r,z)=\psi(z)R(r)

(3)

where $R(r)$ is related only to the range $r$ and satisfies:

\frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}r}\left(r\frac{\mathrm{d}R(r)}{\mathrm{d}r}\right)+k_{r}^{2}R(r)=-\frac{\delta(r)\psi(z_{\mathrm{s}})}{2\pi r\rho(z_{\mathrm{s}})}

(4)

where $k_{r}$ is the horizontal wavenumber. By solving the above formula, we obtain:

R(r)=\frac{\mathrm{i}}{4\rho(z_{\mathrm{s}})}\psi(z_{\mathrm{s}})\mathcal{H}_{0}^{(1)}(k_{r}r)

(5)

where $\mathcal{H}_{0}^{(1)}(\cdot)$ is the Hankel function of the first type and $\psi(z)$ in Eq. (3) satisfies the following modal equation:

\rho(z)\frac{\mathrm{d}}{\mathrm{d}z}\left(\frac{1}{\rho(z)}\frac{\mathrm{d}\psi(z)}{\mathrm{d}z}\right)+k^{2}\psi(z)=k_{r}^{2}\psi(z),\quad k=(1+\mathrm{i}\eta\alpha)\omega/c(z)

(6)

where $k$ is the complex wavenumber, $\alpha$ is the attenuation coefficient in dB $/\lambda$ ( $\lambda$ is the wavelength), and $\eta=(40\pi\log_{10}{\mathrm{e}})^{-1}$ . This is the essential equation to be solved in this paper. When supplemented by boundary conditions, Eq. (6) has a set of solutions $(k_{r,m},\psi_{m}),m=1,2,\dots$ , where $\psi_{m}$ is also called the eigenmode. The eigenmodes of Eq. (6) satisfy orthogonal normalization:

\int_{0}^{H}\frac{{\psi_{m}(z)}{\psi_{n}(z)}}{\rho(z)}\mathrm{d}z=\delta_{mn},\quad m,n=1,2,\dots

(7)

where $H$ is the ocean depth and $\delta$ is the Kronecker delta function. Finally, the fundamental solution of the Helmholtz equation can be written as:

p(r,z)=\frac{\mathrm{i}}{4\rho(z_{\mathrm{s}})}\sum_{m=1}^{\infty}\psi_{m}(z_{\mathrm{s}})\psi_{m}(z)\mathcal{H}_{0}^{(1)}(k_{r,m}r)

(8)

To accurately obtain the sound pressure, it is necessary to synthesize an infinite number of eigenmodes, which is impossible in actual calculations. It is usually more practical to take $M$ physically meaningful eigenmodes to synthesize the sound field. The specific value of $M$ can usually be estimated from the depth of the ocean $H$ , the speed of sound $c$ , and the frequency of the sound source $f$ .

Figure 1: Schematic diagram of the ocean environment with multilayer media.

For ocean environments containing multilayer sediments, $\rho(z)$ and $k(z)$ are usually discontinuous at the interfaces $z=\{h_{\ell}\}_{\ell=1}^{l-1}$ . Considering an intermittent environment, the ocean is divided into $l$ discontinuous layers, as shown by the red dotted line in Figure 1. The environmental parameters are separately defined in the columns as: