Abstract
In this paper, the Chebyshev spectral method is used to solve the normal mode and parabolic equation models of underwater acoustic propagation, and the results of the Chebyshev spectral method and the traditional finite difference method are compared for an ideal fluid waveguide with a constant sound velocity and an ideal fluid waveguide with a deep-sea Munk speed profile. The research shows that, compared with the finite difference method, the Chebyshev spectral method has the advantages of a high computational accuracy and short computational time in underwater acoustic propagation.
1. Introduction
Over 70% of the Earth’s area is covered by the ocean. The ocean is rich in energy, minerals, and biological resources. Countries around the world are vigorously exploring and developing marine resources. Marine protection and development urgently require the detection, identification, positioning, and communication of underwater targets. Due to the characteristics of seawater, electromagnetic waves can be rapidly absorbed, and then they decay in seawater and cannot achieve long-distance transmission. As the main carrier of information transmitted remotely under water, sound waves are of great importance, and it is necessary to deeply understand their propagation features [1, 2].
In mathematical modeling, the law of sound propagation under water is described by a specific set of differential equations. However, because the marine environment is complex and changeable, it is very difficult to obtain an analytical solution of a sound field directly from equations and boundary conditions. Therefore, traditional methods often yield approximate solutions through numerical calculations. To this end, the original differential equations need to be numerically discretized so that the continuous differential equation problem can be converted into computer-solvable discrete algebraic equations. Finally, an approximate solution to the original problem can be obtained. The finite difference method (FDM) has the advantages of being easy to understand and universally applied, so it is still the most common discrete numerical method used in ocean acoustic models.
Common underwater acoustic calculation programs are developed based on the FDM. Normal mode models have been used for many years in underwater acoustics. An early and widely cited paper was proposed by Pekeris [3] who developed a theory for a simple two-layer model of the ocean. Classical normal mode theory cannot address the situation in which the environment changes with distance. The way to address this situation is to segment the distance, and the sound field in each segment is approximately independent of the distance [4]. According to the literature [4–6], the FDM can be used to develop normal mode model for the calculation of an underwater acoustic field, and this approach has become popular among researchers. Based on the parabolic approximation of the original governing equations, previous studies [7–10] used the FDM to develop an approximate calculation method for an underwater acoustic field under wide-angle conditions and developed the numerical calculation software code RAM (Range-dependent Acoustic Model). RAM uses a split-step Padé approximation to solve the parabolic equation; its computational speed and accuracy are greatly improved compared with those of previous parabolic models, so it is widely used. Another study [11, 12] developed a method for an approximate calculation of the underwater acoustic field described by a cylindrical coordinate system using the FDM for a three-dimensional (3D) complex scene and integrated the method into FOR3D. This method has become one of the main methods for studying 3D underwater acoustic propagation.
Although the FDM has been the main method for numerical underwater acoustic field simulations, it still has many shortcomings and deficiencies. For example, the FDM is not flexible enough to handle boundary value problems with complex shapes, and the addition of point and line sources requires special numerical methods; additionally, it is difficult to construct high-precision difference formats, and the convergence and stability of complex numerical computational schemes are difficult to verify mathematically. Notably, the environment must be discretized with a grid that is smaller than the acoustic or elastic wavelength. To meet the rapidly growing needs of various underwater acoustic applications, the development of new discrete methods with high accuracy and rapid calculations has important scientific and application value.
In common scientific and engineering numerical modeling problems, the finite volume method, finite element method, and spectral method are all numerical discrete methods that provide alternatives to the FDM. Among these methods, the finite volume and finite element methods are more suitable for discretizing equations in integral form, and control equations in differential form are more commonly used in the field of underwater acoustics. Therefore, this article explores mainly spectral methods for underwater sound field modeling and calculation applications.
Developed in the 1960s and 1970s, spectral methods are used to numerically solve partial differential equations (PDEs). To date, spectral methods have been widely used in many fields, such as meteorology, physics, and engineering [13–15]. The basic concept of this approach is that the solution of the problem is approximately expanded with a set of smooth basis functions (known as a spectral expansion), thereby transforming the solution problem in the original physical space into an equivalent problem in which the expansion coefficients are determined in the spectral space. In the 1980s, Quartreroni and Canuto systematically studied the theory of the spectral methods [16, 17], concluding that the advantage of the spectral method over other numerical methods in solving linear PDEs is its high accuracy; when solutions of PDEs are smooth enough, the errors of numerical solutions decrease exponentially as the number of discretization nodes increases. Therefore, spectral methods are often regarded as providing “infinite-order” convergence. For periodic boundary conditions, the domain is usually discretized using the Fourier series, while orthogonal polynomials are most frequently used as basis functions for nonperiodic boundary conditions. Among the many types of spectral methods, the Chebyshev spectral method (CSM) with Chebyshev polynomials as the basis functions has been widely used in studies of computational fluid dynamics problems, chemical measurements, and electricity [18–20]. This paper introduces the CSM to solve underwater acoustic sound propagation problems, and the characteristics of the CSM and FDM are compared in terms of the computational accuracy and speed.
2. Underwater Acoustic Propagation Model
In research on underwater acoustic propagation, it is assumed that the seawater density and sound speed are independent of time, and seawater is considered an ideal incompressible fluid. The wave equation in the time-space domain is typically transformed into the Helmholtz equation in the frequency-space domain. The homogeneous Helmholtz equation for a 2D horizontally layered medium with depth and horizontal distance is as follows [21]:where is the sound pressure, and are respectively the density and sound speed of the seawater and are independent of the horizontal distance, and is the angular frequency of the sound source.
Solving the Helmholtz equation in equation (1) directly can lead to various problems, such as those associated with numerical convergence, computation and storage overhead. In practical applications, multiple hypothetical conditions are usually introduced for specific numerical simulation scenarios, and a variety of simplified approximation models, including the normal mode model, parabolic approximation method, and ray method have been developed. Among these models, the normal mode model needs to solve the depth-related equations first to obtain a series of modal solutions and then combines the weighted sum of each mode to obtain the total sound field. The parabolic equation (PE) model is suitable for situations in which the bottom terrain and hydrological parameters change relatively quickly. In this case, the Helmholtz equation can be approximately decomposed into forward waves propagating away from the sound source and backward waves propagating near the sound source, thus ignoring the contribution of reverse acoustic waves to the full wave field; as a result, a parabolic equation model of one-way propagation is obtained [7].
In practical applications, the sound pressure field is usually not used directly; the transmission loss () is generally used, defined aswhere is the sound pressure at a distance of 1 m from the sound source.
2.1. Normal Mode Model
From equation (1), the sound pressure is decomposed into by using the separation of variables method, where satisfies the following equation:where is an unknown constant. Equation (3) is a Sturm-Liouville eigenvalue problem, also known as a modal equation. If supplemented by definitive solution conditions at the sea surface and sea bottom , the modal equation has a series of normal mode solutions . These normal modes are orthogonal to each other and can form a complete set of standard orthogonal bases. The final sound pressure solution obtained by using this property iswhere is the sound source depth and is the first-type Hankel function.
2.2. Parabolic Equation Model
It can be assumed that the solution of equation (1) is given in the form of a cylindrical beam:where is the reference wavenumber; is the envelope function, which represents physical attenuation; and is the Hankel function, which represents spatial dispersion. The Hankel function satisfies the following Bessel differential equation:
Under far-field conditions, that is, , the Hankel function has the following asymptotic form:
By substituting equation (5) into equation (1) and using equations (6) and (7), we obtainwhere is the index of refraction. When the propagation direction is close to the horizontal direction (small horizontal angle), the key paraxial approximation assumption is introduced [7]:
Ignoring the second derivative term, equation (8) becomes
The key to solving the PE is to numerically solve equation (10).
3. CSM for the Underwater Acoustic Propagation Model
The spectral method in this paper is based on using finite-order function expansion and truncation to approximate the original function. The continuous smooth function is expanded according to a set of basis functions in a weighted sum approach. This set of basis functions, also known as trial functions, constitutes a set of orthogonal bases in a continuous function space. If the basis functions are Chebyshev orthogonal polynomials, the CSM is formed. In this case, the trial function is taken as
Any smooth differentiable function defined in the interval can be expanded according to this set of bases as follows:
The expansion coefficient is
Equations (12) and (13) are called the forward Chebyshev transform and inverse Chebyshev transform, respectively.
When numerically evaluating the integral in equation (13), it is necessary to introduce numerical dispersion in the domain . There are many ways to select the integration node . When considering the convenience of the processing boundary conditions and the accuracy of the numerical calculations, the Gauss-Lobatto node , where , is often used. In this case, the expansion coefficient of equation (13) can be approximately calculated as
Additionally, the -order derivative function of can be expanded into
Taking the first derivative as an example, the following relationships are satisfied between the expansion coefficients:
In addition, the spectral coefficients of the product of two known functions can also be represented by the corresponding spectral coefficients:
When applying the CSM, it is necessary to consider the finite truncation of the -term expansion of the infinite term, so the problem of solving an ordinary differential equation is transformed into obtaining a solution for linear algebraic equations with unknowns. How does one start from the control equations and boundary conditions of the physical problem, construct the equations, and solve the equations to obtain the expansion coefficient? Of the available several construction methods, the most common methods are the Tau, Galerkin, collocation, and pseudospectral methods. In recent years, some studies have also aimed to introduce the collocation method and pseudospectral method to acoustic research [22–25]. The authors of these studies also pointed out the advantages and disadvantages of these methods. In the collocation method, linear PDEs are directly solved in physical space so that implementation of the method is easy compared to that of other methods. However, the differentiation matrix derived from the collocation method is full, making computation slow, and ill-conditioning of the differentiation matrix in the collocation method can produce numerical solutions with big errors, especially when a big number of collocation points are used. The pseudospectral method is particularly well suited for solving differential equations with nonlinear terms. The main calculation of this method is still performed in spectral space, but when the method encounters a nonlinear product, the product is first transformed to grid points in physical space, and the result is returned to spectral space to continue the calculation [26–28]. Chowdhury proposed a Rayleigh-Ritz model for the depth eigenproblem of heterogeneous Pekeris waveguides and confirmed the high accuracy of the model through a large number of examples [29]. In this article, we mainly consider the Tau method; this method does not require each trial function to satisfy the boundary conditions, which are directly enforced by the equations of the system. In this approach, the boundary condition is expanded with Eq. (12) and transformed into a linear equation about . Taking the first type of boundary condition as an example, can be expressed as after discretization, and the other boundary conditions can be treated similarly. The constraints of the boundary conditions are reflected in the system of equations used to solve the original problem. In the aforementioned normal mode model and PE of underwater acoustic propagation, the CSM can be introduced to numerically discretize equations (3) and (10) and solve the algebraic eigenvalue problem and the linear equations problem. To this end, a linear transformation is introduced to transform the domain of the original problem from to . Then, equations (3) and (10) become
In the equation, .
3.1. CSM Treatment of the Normal Mode Model
When numerically calculating the normal mode model, we first need to introduce numerical dispersion into equation (3). The traditional method involves performing finite difference discretization in the vertical direction and transforming equation (3) into an algebraic eigenvalue problem to obtain the solution.
The forward Chebyshev transform is used to convert equation (19) from the original physical space to the spectral space; then, by using the -term truncation approximation in conjunction with equations (14) and (18), a linear algebra eigenvalue system is obtained:
Similarly, the boundary condition can be transformed into ; we replace the last two equations in equation (21) to obtain a linear algebraic eigenvalue system in the form of a block:
In the equation, is a square matrix of order , and is a square matrix of order 2. The remaining vectors and divisors of the matrix satisfy the relevant compatibility condition. Equation (22) can be decomposed into the following two equations:
Therefore, it is necessary only to solve the -order algebraic eigenvalue problem in equation (23) and (24). For each set of eigenvalues/eigenvectors (, ), a set of eigensolutions (, ) of equation (21) can be obtained by the inverse Chebyshev transform. Finally, the sound pressure field is obtained by applying equation (4).
3.2. CSM Treatment of the PE
When numerically solving the PE model in equation (10), the traditional method involves using the FDM to discretize and transform the model into a linear equations problem. Here, the operator is used, and the forward Chebyshev transform is used to convert equation (20) from the original physical space to the spectral space. With the -term truncation approximation and equations (14), (17), and (18), we obtain
In the equation, represents the Chebyshev spectral expansion coefficient of the original function , and the operator is the form of in the Chebyshev spectral space. If the semi-implicit format is used for discretization in the direction, then equation (25) can be changed to
Rearranging this equation yieldswhere , , represents the discrete step size in the direction, and represents the expansion coefficient of at grid point .
Similar to the normal mode model, the Tau method is still considered in the PE. The boundary conditions and initial conditions are also transformed to the spectral space by equation (14) and replaced by the last few lines of equation (27).
Equation (27) is a standard -order linear equations problem. After determining the value at , the problem can be solved by a “split-step” calculation to obtain the full field ; then, we use the inverse Chebyshev transform to convert the problem back into the original physical space and obtain the full-field sound pressure solution by equation (5).
A previous study [18] discussed the stability and convergence for the Chebyshev spectral approximation with the Tau method. When a continuity condition and a coercivity condition are assumed for the operator , as in our problems, the stability and convergence can be obtained (Please see the details in [18], chap. 10).
4. Test Case and Analysis of Results
To verify the validity of the CSM in solving the underwater acoustic propagation model, the following tests and analyses are performed through several examples.
Case 1. An ideal fluid waveguide with a constant acoustic velocity
In this simple ocean waveguide model, the density of the seawater is uniform, and the speed of sound is constant. In this case, the sea floor level is considered, the depth is m, the sound source is located at a depth of m, and the sound source frequency is 20 Hz. Pressure release conditions are adopted at the surface and bottom of the sea; that is, . Additionally, the density of the water body is a uniform , and the speed of sound in the water is unchanged at .
According to the wavenumber integration method, the exact analytical solution of the above ideal fluid waveguide sound field isThe vertical wavenumber and horizontal wavenumber are given by the following two formulas:where .
To facilitate the comparison and verification processes, in addition to analytical solution and the CSM, we also tested the Kraken (discrete normal mode model using the FDM) method and the virtual source method (mirror method). Since the PE requires sound waves to propagate in the horizontal direction, this example is not applicable. In Case 1, the receiver is located at a depth of m, and the maximum horizontal distance is 3000 m. For all of the methods, the distance between the discrete grid points in the horizontal direction is taken as 1 m. One thousand equally spaced scattered points are used in the vertical direction for the Kraken. The CSM uses a spectral truncation order of , and the virtual source method uses 500 sound source mirror points. Figure 1 shows the full-field calculated by the analytical solution, virtual source method, Kraken, and CSM. Figure 2 shows the at the receiver’s depth calculated by the four methods. As shown in Figures 1 and 2, the results calculated by the four methods are very similar. The partially enlarged image of Figure 2 reveals that although the results calculated by the four methods are highly similar, there is still a certain error between them in the sound shadow area.
To further investigate the calculation deviations of the various methods, Figure 3 shows the deviations of the other three methods based on the analytical solution. At the receiver’s depth, Kraken produces the largest deviation, followed by the virtual source method, and the region close to the sound source is the singular area of the virtual source method. The CSM produces the smallest deviation from the analytical solution. The error is located mainly in the acoustic shadow areas.
To quantitatively show the accuracy of the CSM and Kraken, we implemented a method described in the literature [29] and calculated the wavenumbers corresponding to the only two modes in this case (equations (29) and (30)). Table 1 shows that the accuracy of the CSM is higher than that of Kraken, regardless of whether the first or second mode is considered.
We implemented different truncation orders and examined the deviations between the calculated results of the CSM and analytical solution. Figure 4 indicates that as the truncation order increases, the deviations in the calculation result decrease rapidly. From to , the deviation decreases by a large margin.
By combining the results of Figures 1–4, we can see that the results of the CSM are correct and reliable. The results also indicate that the CSM needs to use only a small spectral truncation order () to produce satisfactory results. However, the other methods require the use of 1000 discrete grid points in the vertical direction. The CSM has a higher calculation accuracy than Kraken based on the FDM. For the CSM, the calculation error decreases rapidly with additional discrete grid points and then remains at a very low level.
Table 2 shows the computational time of Kraken and the CSM proposed in this paper in Case 1. The results listed in Table 1 are based on Fortran code. In the test, a single thread of a single process is running on a machine node with the Xeon E5-2692 processor of the Tianhe-2 supercomputer. The average times of ten repeated runs are reported. The total run includes two phases: calculating the modes and calculating the sound pressure. From the results of Table 1, the total times of Kraken and the CSM are respectively 43.943 s and 18.738 s, which indicates that the CSM has a faster calculation speed and requires a shorter time. In addition, the CSM takes slightly longer than Kraken to calculate the modes, because the matrix solved by the CSM is dense, while that of Kraken is tridiagonal.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
Case 2. An ideal fluid waveguide with a Munk sound velocity profile
This example is a deep-sea environment example. As shown in Figure 5, the density of the seawater is uniform at . The sound velocity is taken as the Munk sound velocity profile [21], the general form of which iswhere and .
The bottom of the sea is at a depth of 5000 m. We use pressure release conditions at the sea surface and bottom and ignore the effect of the curvature of the Earth on the sound propagation. The sound source is located at a depth of m, and the frequency of the sound source is 50 Hz. The marine hydrological environment in this case can be seen intuitively from Figure 5.
Based on the normal mode model and PE model, we use the CSM and the FDM (the FDM of the normal mode model corresponds to Kraken) for discretization. The receiver is located at a depth of m, and the longest horizontal distance is 100 km. Among these models, the CSM uses a spectral truncation order of , and the distance between discrete points in the vertical direction of the FDM is 1 m; the distance between discrete points in the horizontal direction is 10 m for the two methods. In the Kraken calculation, the range of phase speed is limited to 1500∼1600 m/s, and 103 modes are involved in the calculation of the normal mode model. Since no analytical solution can be given, we compare the relative differences between the 2 models.
Figure 6 shows the full-field calculated using the Kraken, the normal model of the CSM, the PE of the FDM and the PE of the CSM. The results of the four methods are similar, and it is difficult to find any differences between them with the naked eye. To study the differences of various methods more accurately, the curve at the receiver’s depth ( m) calculated by the four methods is given in Figure 7.
Figure 7(a) shows that the results of the two methods based on the normal mode model appear to be uniform (the error between the two is so small that it can be ignored). To quantitatively observe the difference between the FDM and the CSM based on the normal mode model, Table 3 lists the modes (wavenumbers ) calculated by the CSM and Kraken. Since there are 103 modes in Case 2, we list only selected values of in the table for comparison. The wavenumbers calculated by the CSM are very similar to those calculated by Kraken. The difference between the two methods is on the order of .
In Figure 7(b), the initial fields used by the PE model are derived from the normal mode model; the figure illustrates that the results of the two methods based on the PE model are comparable, and the largest difference occurs at a horizontal distance of 60∼80 km; nevertheless, the results are very similar.
The results of Case 2 show that, regardless of whether the CSM is applied to the normal mode model or the PE model, the calculation results are essentially the same as the results of the FDM. Because there is no precise analytical solution for comparison, it is impossible to determine which model is more accurate at slightly different locations, but, in the full field, the results are very similar. However, the significant advantage of the CSM is that it requires fewer orders of discrete algebraic equations to obtain results with an accuracy comparable to that of the FDM results. In this example, the CSM uses only 200 discrete points in the vertical direction, and the results are equivalent to (or even better than) the results obtained with 5000 discrete points in the FDM; thus, the CSM greatly reduces the computational and storage costs. Table 4 shows the computational time of the four methods in Figure 7 (the test platform and method are the same as in Table 2). The computational time of the proposed CSM is usually shorter than that of the FDM in both the normal mode model and PE model, which also confirms our earlier analysis.
The results of the numerical calculation experiments clearly show that the CSM has high accuracy and is an effective method for solving underwater acoustic propagation problems. With the same number of discrete grid points, the CSM can obtain a higher calculation accuracy than other methods. This is especially true when the number of grid points is small. In fact, obtaining the expansion coefficient is equivalent to obtaining an accurate continuous original function, while the FDM requires a larger number of discrete grid points. With approximately the same accuracy, the number of grid points used in the CSM is significantly smaller than that in the FDM. Therefore, the CSM with global approximation ability has obvious advantages in terms of the calculation accuracy and convergence speed.
In the case of the same number of discrete points, the CSM needs to solve the linear equations of the dense matrix vector multiplication problem multiple times when establishing the underwater acoustic propagation model, but the FDM needs only to solve the tridiagonal matrix equations or eigenvalues/eigenvectors problem. Therefore, the calculation time of the CSM is longer than that of the FDM. However, the CSM usually requires fewer discrete points to reach or exceed the accuracy of the FDM. In other words, with the same accuracy, the matrix size required by the CSM is much smaller than that of the FDM. To obtain the spatial distribution of the sound field with the same resolution, the FDM requires more grid points than the CSM; the larger the matrix is, the longer the computational time. However, the computations of the CSM used to solve the underwater acoustic propagation problem require a large number of iterations. The inversion of the dense matrix and requirement of solving equations multiple times result in large memory overhead. Therefore, the CSM and FDM have distinct advantages and disadvantages in terms of the calculation costs. The main disadvantage of the CSM is that the absolute smoothness of the basis function means that it can approximate a smooth original function only. If there are discontinuities or singular points in the distributions of the relevant marine hydrological and environmental parameters, the calculation results will have a large error.
(a)
(b)
(c)
(d)
(a)
(b)
5. Summary and Outlook
In this paper, the CSM is applied to solve the normal mode model and PE model of underwater acoustic propagation, the corresponding mathematical formulas are derived, and the corresponding program codes are developed. The results of the CSM and FDM are compared for an ideal fluid waveguide with a constant acoustic velocity and an ideal fluid waveguide with a deep-sea Munk profile. The results show that the CSM has the advantages of a high computational accuracy and short computational time.
In this article, only the CSM is used for the normal mode and PE models of underwater acoustic propagation. In the future, we plan to extend the spectral method to solve more underwater acoustic calculation models considering cases in which the underwater terrain varies with distance and address the problem of underwater acoustic propagation in more complex underwater environments.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported in part by the National Key Research and Development Program of China (2016YFC1401800), in part by the National Natural Science Foundation of China (61702531, 61972406, and 51709267), and in part by the “Double First-Class” Construction Guidance Project of the National University of Defense Technology (4345161111L).