Next Article in Journal
Acoustic Target Strengths and Swimbladder Morphology of Chub Mackerel Scomber japonicus in the Northwest Pacific Ocean
Previous Article in Journal
An Improved High-Realism Turbulence Simulation of Ocean Scenes in a Maritime Simulator
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Normal Mode Model Based on the Spectral Element Method for Simulating Horizontally Layered Acoustic Waveguides

College of Meteorology and Oceanography, National University of Defense Technology, Changsha 410073, China
*
Author to whom correspondence should be addressed.
J. Mar. Sci. Eng. 2024, 12(9), 1499; https://doi.org/10.3390/jmse12091499
Submission received: 28 July 2024 / Revised: 23 August 2024 / Accepted: 27 August 2024 / Published: 30 August 2024
(This article belongs to the Section Ocean Engineering)

Abstract

:
Acoustic waves are essential tools for guiding underwater activities. For many years, numerical modeling of ocean acoustic propagation has been a major research focus in underwater acoustics. Normal mode theory, one of the earliest and most extensively studied methods in this field, is renowned for its well-established theoretical framework. The core of normal mode theory involves the numerical solution of modal equations. In classical normal mode models, these equations are typically discretized using low-order finite difference methods, which, while broadly applicable, suffer from a limited convergence rate. The spectral element method, widely used in the seismic field, is recognized for its spectral precision and flexibility. In this article, we propose a normal mode model discretized using the spectral element method. The weak form of the modal equation directly satisfies boundary and interface conditions without requiring additional operations. The entire computational domain can be divided into segments of varying number and length, configured according to environmental conditions. The perfectly matched layer technique is employed to simulate acoustic half-space boundary conditions, effectively addressing the high computational costs and numerical instability associated with traditional artificial absorbing layers. Based on these algorithms, we have developed a numerical program (SEM). This research verifies the accuracy of the spectral element model through three different types of numerical experiments.

1. Introduction

Normal mode theory is one of the earliest developed, most widely studied, and theoretically refined methods in underwater acoustics research. It was first applied to underwater acoustics by Pekeris in 1948 [1]. Classical normal mode theory is limited to handling range-independent waveguides. In 1954, Pierce [2] and Miller [3] independently proposed coupled normal mode theory, to address range-dependent waveguides. In 1983, Evans introduced the stepped approximation method, to handle sloping seabeds [4], and in 1986 he developed the numerical model COUPLE, using a new Galerkin variational method to solve the modes [5]. The latest version of this coupled normal mode program, COUPLE07, performs exceptionally well in range-dependent environments, and it is often used as a benchmark solution. Coupled normal mode theory requires consideration of energy exchange between modes, leading to significant computational effort. In the 1960s, Pierce [6] and Milder [7] introduced adiabatic modes for slowly range-dependent waveguides. When environmental fluctuations are gradual, adiabatic modes neglect energy exchange between modes, significantly reducing computational effort.
In recent years, both adiabatic and coupled modes have continued to develop, forming the two main branches of normal mode theory. The key to both approaches lies in solving the local modal equations (or range-independent normal modes). Numerical methods, such as finite element, finite difference, and spectral methods, are typically used to convert the continuous modal equations into algebraic equations that can be solved on a computer.
The well-established models include KRAKEN [8], which utilizes finite difference discretization, and COUPLE, which employs Galerkin variational discretization. The finite difference method offers several advantages: it is simple, intuitive, easy to understand, performs well in various complex environments, and has strong versatility and numerical stability. Additionally, it features straightforward grid setup and quick generation. However, low-order finite difference schemes have slow convergence rates, while high-order schemes can be inconvenient to implement. On the other hand, the Galerkin method provides high accuracy and fast convergence, but it requires pre-constructing basis functions that satisfy boundary and interface conditions, making it challenging to apply in scenarios with complex boundary conditions.
This paper proposes a normal mode model that is discretized using the spectral element method. The spectral element method is an approach that combines the advantages of spectral methods and finite element methods [9,10,11]. Initially developed in computational fluid dynamics, it was first introduced by Patera in 1984, employing Chebyshev polynomial expansions within elements, to obtain numerical solutions to the Navier–Stokes equations [12]. Essentially, it can be viewed as a hybrid of spectral methods and finite element methods, offering the precision of spectral methods and the flexibility of finite element methods [13].
In 1994, Priolo and Carcione successfully applied this method to simulate the propagation of acoustic and elastic waves [14]. In 1989, Ho and Maday introduced Lagrange interpolation basis functions into the spectral element method, combining them with Gauss–Legendre–Lobatto quadrature. By using the integration points as interpolation grids, they developed the Lagrange spectral element method [15]. Seriani and colleagues applied the spectral element method to wave propagation in range-independent media and three-dimensional acoustic modeling in 1995 and 1998 [16]. In the early 21st century, many researchers conducted detailed theoretical analyses of the spectral element method, proving its accuracy and convergence [17,18,19,20]. In 2011, Peter et al. introduced its application in seismic wave propagation [21]. In 2012, Cristini and Komatitsch used the spectral element method to simulate two-dimensional acoustic fields, comparing the results with those from the SCOOTER [22] program and demonstrating its potential in underwater acoustics [23]. Currently, research and applications of the spectral element method, both domestically and internationally, are primarily focused on the seismic field [24,25]. Given the similarities between seismology and ocean acoustics, some have attempted to use SPECFEM for time-domain problems in ocean acoustics. SPECFEM2D and SPECFEM3D (https://github.com/) are open-source software packages based on the spectral element method for simulating seismic wave propagation in the time domain. They have been widely used in seismology, petroleum exploration, geophysical research, and other industries, becoming essential tools in these fields. In recent years, Shi proposed a domain decomposition method based on the spectral element method for elastic wave computation in the frequency domain, confirming its spectral accuracy [26,27]. However, the spectral element method has not yet been widely applied to frequency-domain problems in ocean acoustics.
For this paper, we employed the spectral element method to discretize the modal equations of the normal mode model, and we developed the SEM model for simulating range-independent ocean acoustic waveguides. SEM was used to perform numerical simulations on pseudolinear-speed waveguides, waveguides with two-layer media, and Munk waveguides. The experimental results were compared with those from classic underwater acoustic simulation software, such as KRAKEN and COUPLE, demonstrating the accuracy of SEM. SEM represents the initial exploration of applying the spectral element method to frequency-domain ocean acoustics. Moreover, SEM integrates the perfectly matched layer (PML) technique, to effectively handle acoustic half-space boundaries. The simulation results also validate the effectiveness of the SEM.

2. Theory and Model

2.1. Normal Mode Model

For the ocean environment shown in Figure 1, the properties of seawater remain constant in the horizontal direction. The waveguide environment in the vertical direction consists of multiple layers of media, with the acoustic parameters in each layer defined as follows:
c ( z ) = c l ( z ) , h l 1 z h l , l = 1 , , L c , z H
ρ ( z ) = ρ l ( z ) , h l 1 z h l , l = 1 , , L ρ , z H
α ( z ) = α l ( z ) , h l 1 z h l , l = 1 , , L α , z H
where c ( z ) is the sound speed of the media, ρ ( z ) is the density, and α ( z ) is the attenuation in d B / λ . First, we consider a simple example where the acoustic properties of the media are range-independent. We start from the Helmholtz equation for acoustic pressure:
1 r r r p r + ρ ( z ) z 1 ρ ( z ) p z + k 2 ( z ) p = δ ( r ) δ z z s 2 π r
where k ( z ) = 2 π f / c ( z ) . The range independence allows us to use the separation-of-variables technique to decompose the acoustic pressure into a depth function and a range function: p ( r , z ) = m = 0 M ψ m ( z ) R m ( r ) , where the depth function ψ m ( z ) satisfies the modal equation:
d d z 1 ρ ( z ) d ψ m ( z ) d z + k 2 ( z ) ψ m ( z ) ρ ( z ) = k r , m 2 ψ m ( z ) ρ ( z )
The modal equation, Equation (5), is a classical Sturm–Liouville eigenvalue problem. Solving this eigenvalue problem directly yields multiple pairs of eigenvalues and eigenvectors, where the square root of the eigenvalue is referred to as the horizontal wavenumber, and the eigenvector is the mode function. The modes satisfy the orthogonal normalization condition:
0 H ψ m 1 ( z ) ψ m 2 ( z ) ρ ( z ) d z = 1 , m 1 = m 2 0 , m 1 m 2
The range function R m ( r ) after the separation of variables can be expressed using the Hankel function of the first kind:
R m ( r ) = i 4 ρ ( z s ) ψ m ( z s ) H 0 ( 1 ) ( k r , m , r )
Thus, the acoustic pressure solution of the classical normal mode theory is obtained:
p ( r , z ) = i 4 ρ z s m = 1 M ψ n z s ψ m ( z ) H 0 ( 1 ) k r , m , r

2.2. Perfectly Matched Layer Technique

In ocean acoustic problems, modeling the waveguide with a homogeneous half-space as the lower boundary is a classic approach. However, normal mode theory employs the residue theorem to convert the continuous integral problem into a discrete modal solution. When the lower boundary is an acoustic half-space, the modal equation becomes a singular Sturm–Liouville problem, which is challenging to solve with conventional methods. Existing normal mode models address the acoustic half-space condition by using artificial absorbing layers to mitigate the clutter reflected from the seabed (as in COUPLE), simulating the attenuation of acoustic waves in an infinitely deep half-space. However, for low-frequency acoustic field simulations, the thickness of the artificial absorbing layer can greatly exceed the depth of the waveguide, resulting in high computational costs.
The PML technique effectively addresses these issues. Introduced by Berenger in 1994 [28], PML was designed to overcome the significant computational slowdowns associated with thick artificial absorbing layers. Since its inception, PML has been widely used in acoustics, electromagnetics, geophysical fluid dynamics, quantum mechanics, and other fields [29,30]. A well-designed PML can absorb reflected waves and simulate a half-space with a thickness of only a few wavelengths [31,32].
Traditional artificial absorbing layers use medium attenuation coefficients to suppress reflected waves. In contrast, the key to the PML technique lies in complex coordinate transformation, which attenuates reflected acoustic waves through coordinate damping:
z ^ = z + i 0 z V ( x ) d x
V ( z ) = 0 , 0 z < H β z H D H 2 , H z D
where ( D H ) is the thickness of the perfectly matched layer and β represents the PML absorption coefficient, controlling the absorption effectiveness of the PML. With an appropriate absorption coefficient setting, the wave can be attenuated to zero within a thickness of two wavelengths. When the PML is present, the modal equation with the coordinate transformation becomes [33]:
1 σ ( z ) d d z 1 σ ( z ) 1 ρ ( z ) d ψ ( z ) d z + k 2 ( z ) ψ ( z ) ρ ( z ) = k r 2 ψ ( z ) ρ ( z )
σ ( z ) : = d z ^ d z = 1 + i V ( z )

3. Spectral Element Method and Numerical Discretization

3.1. Spectral Element Method

The spectral element method combines the high accuracy and fast convergence of spectral methods with the flexibility of finite element methods, making it easier to handle complex environments. In practical numerical discretization, the spectral element method first discretizes the solution domain into several elements, as does the finite element method, and it then applies spectral discretization within each element. Finally, the discretized linear equations within each element are solved simultaneously. In the spectral discretization process within each element, basis functions need to be introduced. A set of N th -order Lagrange polynomials { n ( x ) } n = 0 N associated with the given nodes x i on the domain are used as basis functions. The Lagrange interpolation polynomial corresponding to the node x i is denoted as i ( x ) , which is defined by all the nodes:
n ( x ) : = i = 0 , i n N x x n x i x n , i = 0 , 1 , 2 , , N
The function u ( x ) can be expanded as a linear combination of Lagrange basis functions:
u N ( x ) = j = 0 N u N ( x j ) j ( x )
The coefficients of the expansion terms are related to the selection of the grid points. The Lagrange polynomials are equivalent to the Dirac function δ :
i ( x j ) : = δ i j = 1 , i = j 0 , i j , i = 0 , 1 , 2 , , N
That is, n ( x ) has a value only at the node x n , and it is zero at all other nodes. Equation (5) is strictly satisfied at the chosen grid points. For the convenience of subsequent operations, such as integration, differentiation, and multiplication, this paper uses nodes x 0 = 1 , x N = 1 , x k ( 0 < k < N ) as Gauss–Legendre–Lobatto (GLL) grid points. These points are the projections of the stationary points of the N th -order Legendre polynomial T N ( x ) onto the domain of the equation to be solved. GLL grid points are characterized by a dense distribution at the ends and a sparse distribution in the middle:
T n ( x ) = 1 2 n n ! d n d x n [ ( x 2 1 ) n ]
In the grid space, the product f ( x ) of two functions u ( x ) and v ( x ) can be directly expressed as:
f ( x j ) : = u ( x j ) v ( x j )
To facilitate the subsequent representation of the product of two functions, we define a matrix C f :
C k j = f ( x j ) k = j ; 0 , k j .
At this point, C f can be represented as C u C v . C f is composed of the values of the corresponding function at the respective nodes; hence, the order of the C f matrix is the same as the number of nodes. Using the function values at the grid points, the integral can be obtained using the following Gauss–Legendre–Lobatto quadrature formula:
1 1 f ( x ) d x j = 0 N f ( x j ) w j , w j : = 2 N ( N + 1 ) T N ( x j ) 2
The first-order differentiation matrix D N corresponding to N th nodes can be expressed as:
D k j = T N ( x k ) ( x k x j ) T N ( x j ) , k j ; N ( N + 1 ) 4 , k = j = 0 ; N ( N + 1 ) 4 , k = j = N ; 0 , others .

3.2. Numerical Discretization

Next, we introduce the spectral element discretization of the modal equation. First, we define a test function w ( x ) . By multiplying both sides of the modal equation, Equation (5), by the test function and taking the inner product, we obtain the weak form:
0 H d d x 1 ρ ( x ) d ψ ( x ) d x w ( x ) + κ 2 ( x ) ρ ( x ) ψ ( x ) w ( x ) d x = 0 H κ r 2 ρ ( x ) ψ ( x ) w ( x ) d x
We performing integration by parts on the first term of the equation:
0 H d d x 1 ρ ( x ) d ψ ( x ) d x w ( x ) d x = 0 H d d x w ( x ) ρ ( x ) d ψ ( x ) d x d x 0 H 1 ρ ( x ) d w ( x ) d x d ψ ( x ) d x d x
Given the boundary conditions p = 0 and ψ = 0 , the weak form of the modal equation is finally obtained as
0 H k 2 ( z ) ρ ( z ) ψ ( z ) w ( z ) 1 ρ ( z ) d w ( z ) d z d ψ ( z ) d z d z = 0 H k r 2 ρ ( z ) ψ ( z ) w ( z ) d z
From a numerical-solution perspective, the weak form of the modal equation inherently satisfies the boundary conditions. Thus, using the weak form not only effectively simulates wave propagation in elastic media but also enhances accuracy for simulations involving inelastic media, surface topography, and interface wave propagation. In an actual vertically inhomogeneous ocean medium, the region is divided into E computational elements [ h e 1 , h e ] . When the environment is inherently layered, a practical approach is to place an element above and below each layer interface. For each layer of the medium, computational elements can be tailored based on environmental complexity. Within each computational element, suitable collocation points are selected to form local Lagrange interpolation basis functions, which are then combined to create global basis functions.
In Figure 2. above, end = 1 + e = 1 E N e . Note that adjacent Element e and Element e + 1 share collocation points and interpolation basis functions at their interface. According to the continuity condition, the local degrees of freedom N e e and 0 e + 1 represent the same global degree of freedom. Therefore, the global basis functions can be expressed as follows: By the continuous condition that the local degree of freedom corresponds to, N e e and 0 e + 1 refer to the same global degree of freedom. So, the global basis L set is 0 1 , 1 1 , 2 1 , , N 1 1 1 , N 1 1 + 0 2 , 1 2 , , N E E , where it is defined as follows:
i e ( x ) = i e ( x ) , x Element e , 0 , others .
By taking the test function w ( x ) in Equation (23) as any basis function k e ( x ) in the L space, we note that k e ( x ) is only defined within Element e; thus, we obtain:
h e 1 h e k 2 ( z ) ρ ( z ) ψ ( z ) k e ( z ) 1 ρ ( z ) d k e ( z ) d z d ψ ( z ) d z d z = h e 1 h e k r 2 ρ ( z ) ψ ( z ) k e ( z ) d z
Since GLL nodes are defined on the interval [ 1 , 1 ] , the following mapping relationship is used to map the solution domain from [ h e 1 , h e ] to the interval [ 1 , 1 ] :
x = 2 z h e h e 1 h e + h e 1 h e h e 1 , J = d x d z = 2 h e h e 1 = 2 Δ H e
We define ψ j e = ψ e ( x j ) to represent the value of the modal function at the jth collocation point in Element e, and we expand ψ ( z ) using Equation (14):
j = 0 N e ψ j e 1 1 k e ( x ) 2 ρ ( x ) k e ( x ) j e ( x ) 4 Δ H e 2 1 ρ ( x ) d k e ( x ) d x d j e ( x ) d x d x = j = 0 N e ψ j e 1 1 k r 2 ρ ( x ) k e ( x ) j e ( x ) d x
Using the Gaussian quadrature formula, Equation (19), to convert the integrals in Equation (27) into a summation form, and noting that k e ( x j ) = 0 when k j , the first term on the left-hand side and the first term on the right-hand side of the equation have values only when j = k :
ψ k e k e ( x k ) 2 ρ ( x k ) w k e j = 0 N e 4 ψ j e Δ H e 2 i = 0 N e w i e ρ ( x i ) d k e ( x ) d x d j e ( x ) d x x = x i = ψ k e k r 2 ρ ( x k ) w k e
Ψ e is the column vector consisting of ψ j e j = 0 N e . The ( N e + 1 ) equations in Element e can be represented in matrix form as:
A e Ψ e + B e Ψ e = Q e Ψ e
The matrix relationship and the expressions for the elements of each matrix are shown below:
A e = C w e C k 2 e C 1 / ρ e , A k j e = k e ( x k ) 2 ρ ( x k ) w k e , k = j 0 , k j
B e = 4 Δ H e 2 D N e C w e C 1 / ρ e D N e T , B k j e = 4 Δ H e 2 i = 0 N e w i e ρ ( x i ) d k e ( x ) d x d j e ( x ) d x x = x i
Q e = k r 2 C w e C 1 / ρ e , Q k j e = k r 2 ρ ( x k ) w k e , k = j 0 , k j
When the lower boundary is an acoustic half-space, an additional computational element for the PML is required. According to Equation (11), the control matrices in the PML layer can be obtained by writing C 1 / ρ E + 1 in A and Q as C 1 / ρ E + 1 C σ E + 1 and by writing C 1 / ρ E + 1 in B as C 1 / ρ E + 1 C 1 / σ E + 1 . The global matrices A , B , and Q are obtained by diagonally assembling the submatrices of all the computational elements:
A = A 1 0 0 0 0 A 2 0 0 0 0 0 0 0 0 A E , B = B 1 0 0 0 0 B 2 0 0 0 0 0 0 0 0 B E Q = Q 1 0 0 0 0 Q 2 0 0 0 0 0 0 0 0 Q E
For the adjacent Element e and Element e + 1 , the interface conditions should be satisfied:
ψ N e e = ψ 0 e + 1
1 ρ N e e d ψ N e e d z = 1 ρ 0 e + 1 d ψ 0 e + 1 d z
After dividing the computational elements, simply substituting the boundary conditions in Equation (23) is not strictly valid:
0 H d d x w ( x ) ρ ( x ) d ψ ( x ) d x d x = w ( x ) ρ ( x ) d ψ ( x ) d x | 0 H + e = 1 E 1 w ( h e ) 1 ρ N e e d ψ N e e d z 1 ρ 0 e + 1 d ψ 0 e + 1 d z
By substituting the boundary conditions and the particle velocity continuity condition Equation (35) into Equation (36), Equation (36) can still be simplified to 0, which means Equation (23) holds strictly. Additionally, we must ensure that the governing equations at the interface satisfy the pressure continuity condition at the interface:
M N e N e e ψ N e e + j = 0 N e K N e j ψ j e = B N e N e e ψ N e e
M 00 e + 1 ψ 0 e + 1 + j = 0 N e + 1 K 0 j ψ j e + 1 = B 00 e + 1 ψ 0 e + 1
j = 0 N e 1 K N e j ψ j e + K N e N e e + K 00 e + 1 ψ N e e + j = 1 N e + 1 K 0 j ψ j e + 1 = M N e N e e + M 00 e + 1 ψ N e e + B N e N e e + B 00 e + 1 ψ N e e
By using Equation (39) to replace Equation (37) and Equation (38), the resulting equations naturally satisfy the pressure continuity condition. Correspondingly, we replace the two rows of equations at the interfaces in A , B , and Q with Equation (39). After processing the E − 1 interfaces, the final governing equations are obtained:
A Ψ + B Ψ = Q Ψ
It is noted that at this point Ψ = ψ ( z 0 ) , ψ ( z 1 ) , ψ ( z 2 ) , , ψ ( z end 1 ) , ψ ( z end ) T , representing the values of the modal functions at all collocation points. Considering that the boundary conditions provide the function values at the first and last nodes, solving the ( end 2 ) order generalized matrix eigenvalue problem yields the corresponding ( end 2 ) sets of eigenvalues and eigenvectors k r , m 2 , ψ m m = 0 end 2 . We select the required modes and normalize; then, the final acoustic pressure can be calculated, using Equation (8).

4. Numerical Simulation

4.1. Pseudolinear-Speed Waveguide

The pseudolinear-speed waveguide is an analytical example, where the acoustic field has an analytical solution. The distribution of sound speed with depth in the pseudolinear-speed waveguide satisfies the following equation:
c ( z ) = 1 a z + b
In the example, the parameter was set to a = 2.78 × 10 10 s 2 / m 3 , b = 4.2 × 10 7 s 2 / m 2 . The ocean depth was H = 100 m, the density was set to ρ = 1 g/cm3, and the seabed was modeled as a pressure-free bottom. The source depth was z s = 36 m, with a frequency of f = 50 Hz. In this example, SEM equidistantly divided the computational domain into 8 elements, each using 7th-order GLL nodes.
Figure 3a–c show the acoustic fields of the pseudolinear-speed waveguide with a pressure-free bottom, as calculated by the analytical solution, KRAKEN, and SEM. The number of discrete points used was 100 for KRAKEN and 48 for SEM. The global error of SEM was 0.0064 dB, while KRAKEN had an error of 0.0168 dB. The transmission loss (TL) curves at a depth of z = 40 m versus range, shown in Figure 3d, indicate that the results from SEM and KRAKEN at this depth are consistent with the analytical solution. Table 1 lists the horizontal wavenumbers of the six modes used to synthesize the acoustic field in both KRAKEN and SEM. The value of k r represents the horizontal wavenumber and corresponds to the eigenvalues obtained from the eigenvalue problem in the discretized modal equations. The first seven significant digits of the horizontal wavenumbers calculated by SEM match those obtained by KRAKEN, demonstrating the reliability of the modes calculated by SEM. Comparison of the transmission loss curves and average error results shows that the acoustic fields simulated by SEM and KRAKEN closely match the analytical solution, confirming the high accuracy of SEM and validating the correctness of the SEM model.

4.2. A Waveguide with Two-Layer Media

To assess the accuracy of SEM in layered environments with significant vertical variations in sound speed, we examined a waveguide with two-layer media for this subsection. The marine environment was modeled with two layers: the upper layer, extending to 800 m depth, exhibited significant variations in sound speed with a density set to ρ = 1 g/cm3. Below 800 m, there was a sedimentary layer. The sound speed distribution is shown in Figure 4a. In this example, SEM divided the computational domain into 10 elements, each using 15th-order GLL nodes.
Figure 5 displays the mode shapes for the two-layer waveguide, as calculated by KRAKENC (blue lines) and SEM (red dashed lines) at 33 Hz. KRAKENC, an updated version of KRAKEN [8], uses root-finding methods to determine the eigenvalues of the singular Sturm–Liouville problem. Although KRAKENC offers improved accuracy over KRAKEN, it also comes with increased runtime. In this example, the modal functions obtained from SEM and KRAKENC were nearly identical, with both methods showing high agreement in the detailed variations of all modes. Figure 6a,c present the acoustic fields for the two-layer waveguide with a pressure-free bottom, as calculated by KRAKENC, COUPLE, and SEM. In the first layer, the simulation results are highly consistent, accurately capturing both the near-field convergence zone and the far-field reflected waves. However, in the sediment layer, there were subtle differences between the three models’ calculations. Specifically, in mode 23, the models diverged below 800 m. Additionally, the TL curve shown in Figure 6d at a depth of 100 m nearly overlaps perfectly. In this case, COUPLE remained consistent with SEM throughout, while KRAKENC showed a difference only in the extremum calculation at r = 1000 m.

4.3. Munk Waveguide

The Munk waveguide represents a classic deep-sea example. In this case, the ocean depth was H = 5000 m, and the density was set to ρ = 1 g/cm3. The waveguide environment is depicted in Figure 4b. Below 5000 m, the sound speed abruptly changed to 1600 m/s, while density and attenuation remained unchanged. For this example, SEM divided the computational domain equidistantly into 20 elements, each utilizing 30th-order GLL nodes.
Figure 7a–c show the acoustic fields of the Munk waveguide with an acoustic half-space, as calculated by KRAKENC, COUPLE, and SEM. COUPLE uses artificial absorbing layers to simulate the acoustic half-space, and SEM utilizes the PML technique ( β = 5 ( 6 H f / c ) 0.25 ) to model the acoustic half-space. In terms of far-field results, all three models accurately computed the propagating modes. However, for near-field leaky modes, COUPLE performed less effectively compared to KRAKENC and SEM, as COUPLE failed to compute the lateral waves within the near-field dashed box. The results from KRAKENC and SEM are more consistent and better align with the physical principles of acoustic field propagation. As shown in Figure 7d, the TL results from SEM are highly consistent with those from KRAKENC and COUPLE, both in near-field and far-field conditions. The nearly identical TL curves suggest that SEM can simulate long-range deep-sea waveguides with performance comparable to KRAKENC and COUPLE. Figure 8 displays the horizontal wavenumbers for 329 modes calculated by SEM and KRAKENC. For modes with small real parts of the horizontal wavenumbers, there are slight differences in the imaginary parts between SEM and KRAKENC. However, for modes with large real parts, the results are essentially identical. This suggests that while the near-field leaky modes exhibit minor differences between the two models, the propagative modes in the far field are very similar. This agreement is also observed in the modal functions. Modes are indexed based on the magnitude of the real part of the horizontal wavenumber. The first three modes in Figure 9 show that the primary mode shapes calculated by KRAKENC closely match those from SEM. However, the fourth subfigure, which corresponds to a higher-order mode with a smaller real part and a significant imaginary part, shows some distortion. This is evident from the deviations in the peak values of the mode shapes calculated by KRAKENC and SEM below 3000 m.

5. Discussion

In the numerical simulations described, we present the values of k r , compare the TL at specific depths, and examine the mode shapes. These examples demonstrate the accuracy of SEM across various ocean environments: simple ocean environments, layered environments with significant sound speed changes, and deep-sea conditions. In the analytical example of the pseudolinear-speed waveguide, where the acoustic properties vary relatively simply, SEM achieved high accuracy with fewer collocation points. For the waveguide with two-layer media, where the sound speed varies dramatically within 800 m, SEM’s results are highly consistent with those from KRAKEN and COUPLE. As shown in Figure 2, SEM divided the entire domain into multiple computational elements, avoiding the issue of the sparse collocation points that can occur with a single-element setup.
For specific example settings, environmental parameters can be adjusted by increasing the number of computational elements or collocation points in regions with significant changes. For the waveguide with two-layer media discussed in this paper, we can use either 100 equal-length computational elements with 3rd-order grid points or 10 equally spaced elements with 15th-order grid points. Both approaches allow for accurate computation of the acoustic field. Refining the model by increasing the number of elements or collocation points in regions with significant acoustic parameter variations enables more detailed and accurate acoustic field information tailored to the specific needs of the simulation. In the classic deep-sea example of the Munk waveguide, the PML technique effectively manages the acoustic half-space. COUPLE uses a 3000 m thick artificial absorbing layer (50 wavelengths), whereas SEM, employing the PML technique, achieves equivalent functionality with only a 45 m thick layer (0.75 wavelength). For deep-sea examples, the computational cost of PML is minimal. Integrating PML technique for handling the acoustic half-space significantly conserves computational resources, allowing more focus on areas of interest. Under acoustic half-space conditions, both KRAKEN and COUPLE encounter issues with root-missing when using root-finding algorithms to compute eigenvalues. Therefore, KRAKENC, which offers better numerical stability, was used in this example. A detailed comparison of horizontal wavenumbers and modal functions calculated by SEM and KRAKENC confirms the accuracy of the normal modes calculated by SEM, especially after incorporating the PML technique.

6. Conclusions

This paper introduced an SEM model designed to address the challenge of acoustic propagation in vertically stratified ocean environments, incorporating the PML technique. Our numerical simulations demonstrated that SEM accurately calculates acoustic fields, while the PML technique effectively manages acoustic half-space conditions. The flexibility in dividing computational elements and adjusting the order of grid points enhances the model’s adaptability. Although this study focused on range-independent cases and used relatively simple numerical examples, SEM marks a significant advancement in discretizing and solving frequency-domain acoustic models using the spectral element method. The simulation results confirm the feasibility of SEM for frequency-domain underwater acoustics and highlight the potential for further theoretical exploration. Future research will aim to develop 2D and 3D frequency-domain spectral element models, incorporating more realistic boundary conditions such as elastic boundaries, and integrating these models with experimental data for further investigation.

Author Contributions

Conceptualization, Y.W., H.T. and G.X.; methodology, Y.Z. and H.T.; software, Y.Z.; validation, Y.Z. and D.G.; formal analysis, Y.Z.; writing, Y.Z. and H.T.; funding acquisition, Y.W. and G.X. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported in part by the Hu’nan Provincial Natural Science Foundation (Grant No. 2024JJ5409).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Pekeris, C.L. Theory of propagation of explosive sound in shallow water. Geol. Soc. Am. Mem. 1948, 27, 1–117. [Google Scholar] [CrossRef]
  2. Pierce, J.R. Coupling of Modes of Propagation. J. Appl. Phys. 1954, 25, 179–183. [Google Scholar] [CrossRef]
  3. Miller, S.E. Coupled wave theory and waveguide applications. Bell Syst. Tech. J. 1954, 33, 661–719. [Google Scholar] [CrossRef]
  4. Evans, R.B. A coupled mode solution for acoustic propagation in a waveguide with stepwise depth variations of a penetrable bottom. J. Acoust. Soc. Am. 1983, 74, 188–195. [Google Scholar] [CrossRef]
  5. Evans, R.B. The decoupling of stepwise coupled modes. J. Acoust. Soc. Am. 1986, 80, 1414–1418. [Google Scholar] [CrossRef]
  6. Pierce, A.D. Extension of the Method of Normal Modes to Sound Propagation in an Almost-Stratified Medium. J. Acoust. Soc. Am. 1965, 37, 19–27. [Google Scholar] [CrossRef]
  7. Milder, D.M. Ray and Wave Invariants for SOFAR Channel Propagation. J. Acoust. Soc. Am. 1969, 46, 1259–1263. [Google Scholar] [CrossRef]
  8. Porter, M.B. The Kraken Normal Mode Program; SACLANT Undersea Research Centre: La Spezia, Italy, 2001. [Google Scholar]
  9. Tu, H.; Wang, Y.; Zhang, Y.; Wang, X.; Liu, W. A spectrally discretized wide-angle parabolic equation model for simulating acoustic propagation in laterally inhomogeneous oceans. J. Acoust. Soc. Am. 2023, 153, 3334–3349. [Google Scholar] [CrossRef]
  10. Tu, H.; Wang, Y.; Yang, C.; Liu, W.; Wang, X. A Chebyshev–Tau spectral method for coupled modes of underwater sound propagation in range-dependent ocean environments. Phys. Fluids 2023, 35, 037113. [Google Scholar] [CrossRef]
  11. Tu, H.; Wang, Y.; Liu, W.; Ma, S.; Wang, X. A spectral method for the depth-separated solution of a wavenumber integration model for horizontally stratified fluid acoustic waveguides. Phys. Fluids 2023, 35, 057127. [Google Scholar] [CrossRef]
  12. Patera, A.T. A spectral element method for fluid dynamics: Laminar flow in a channel expansion. J. Comput. Phys. 1984, 54, 468–488. [Google Scholar] [CrossRef]
  13. Seriani, G.; Carcione, J.M.; Priolo, E. Numerical simulation of interface waves by high-order spectral modeling techniques. J. Acoust. Soc. Am. 1992, 92, 2456. [Google Scholar] [CrossRef]
  14. Weijun, L. A Chebyshev spectral element method for elastic wave modeling. Acta Acust. 2007, 42, 525–533. [Google Scholar] [CrossRef]
  15. Lee-Wing, H.; Maday, Y.; Patera, A.T.; Rønquist, E.M. A high-order Lagrangian-decoupling method for the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Eng. 1990, 80, 65–90. [Google Scholar] [CrossRef]
  16. Seriani, G. 3-D large-scale wave propagation modeling by spectral element method on Cray T3E multiprocessor. Comput. Methods Appl. Mech. Eng. 1998, 164, 235–247. [Google Scholar] [CrossRef]
  17. Deville, M.O.; Fischer, P.F.; Mund, E.H.; Gartling, D.K. High-Order Methods for Incompressible Fluid Flow. Appl. Mech. Rev. 2003, 56, B43. [Google Scholar] [CrossRef]
  18. De Basabe, J.D.; Sen, M.K. Grid dispersion and stability criteria of some common finite-element methods for acoustic and elastic wave equations. Geophysics 2007, 72, T81–T95. [Google Scholar] [CrossRef]
  19. Cohen, G.C.; Liu, Q.H. Higher-Order Numerical Methods for Transient Wave Equations. J. Acoust. Soc. Am. 2003, 114, 21. [Google Scholar] [CrossRef]
  20. Melvin, T.; Staniforth, A.; Thuburn, J. Dispersion analysis of the spectral element method. Q. J. R. Meteorol. Soc. 2012, 138, 1934–1947. [Google Scholar] [CrossRef]
  21. Peter, D.; Komatitsch, D.; Luo, Y.; Martin, R.; Le Goff, N.; Casarotti, E.; Le Loher, P.; Magnoni, F.; Liu, Q.; Blitz, C.; et al. Forward and adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes. Geophys. J. Int. 2011, 186, 721–739. [Google Scholar] [CrossRef]
  22. Porter, M.B. SCOOTER: A Finite Element FFP Code. 2001. Available online: https://oalib-acoustics.org/models-and-software/acoustics-toolbox/ (accessed on 26 August 2024).
  23. Cristini, P.; Komatitsch, D. Some illustrative examples of the use of a spectral-element method in ocean acoustics. J. Acoust. Soc. Am. 2012, 131, EL229–EL235. [Google Scholar] [CrossRef] [PubMed]
  24. Komatitsch, D.; Tromp, J. Spectral-element simulations of global seismic wave propagation—I. Validation. Geophys. J. Int. 2002, 149, 390–412. [Google Scholar] [CrossRef]
  25. Tongkui, W.; Ruihua, L.; Xiaofan, L.; Meigen, Z.; Guihua, L. Numerical spectra-l element modeling for seismic wave propagation in transversely isotropic medium. Prog. Geophys. 2007, 22, 778–784. [Google Scholar]
  26. Komatitsch, D.; Tromp, J.; Shi, L.; Zhou, Y.; Wang, J.; Zhuang, M.; Liu, N.; Liu, H. Spectral element method for elastic and acoustic waves in frequency domain. J. Comput. Phys. 2016, 327, 19–38. [Google Scholar] [CrossRef]
  27. Shi, L.; Zhuang, M.; Zhou, Y.; Liu, N.; Liu, Q. Domain decomposition based on the spectral element method for frequency-domain computational elastodynamics. Sci. China Earth Sci. 2021, 64, 388–403. [Google Scholar] [CrossRef]
  28. Berenger, J.P. A perfectly matched layer for the absorption of electromagnetic waves. J. Comput. Phys. 1994, 114, 185–200. [Google Scholar] [CrossRef]
  29. Chew, W.C.; Weedon, W.H. A 3D perfectly matched layer for the absorption of electromagnetic waves. Microw. Opt. Technol. Lett. 1994, 7, 599–604. [Google Scholar] [CrossRef]
  30. Lu, Y.Y.; Zhu, J. Perfectly matched layer for acoustic waveguide modeling — benchmark calculations and perturbation analysis. Comput. Model. Eng. Sci. 2007, 22, 235–248. [Google Scholar] [CrossRef]
  31. Singer, I.; Turkel, E. A perfectly matched layer for the Helmholtz equation in a semi-infinite strip. J. Comput. Phys. 2004, 201, 439–465. [Google Scholar] [CrossRef]
  32. Rabinovich1, D.; Givoli1, D.; Bécache, E. Comparison of high-order absorbing boundary conditions and perfectly matched layers in the frequency domain. Int. J. Numer. Methods Biomed. Eng. 2010, 26, 1351–1369. [Google Scholar] [CrossRef]
  33. He, T.; Mo, S.; Guo, W.; Fang, E. Modeling propagation in shallow water with the range-dependent sea surfaces and fluid seabeds using the equivalent source method. J. Acoust. Soc. Am. 2021, 149, 997–1011. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Schematic diagram of vertically inhomogeneous oceans.
Figure 1. Schematic diagram of vertically inhomogeneous oceans.
Jmse 12 01499 g001
Figure 2. Grid configuration points and basis functions.
Figure 2. Grid configuration points and basis functions.
Jmse 12 01499 g002
Figure 3. Acoustic fields of the pseudolinear-speed waveguide with a pressure-free bottom calculated by analytical solution (a), KRAKEN (b), and SEM (c); TL curves at a depth of z = 40 m versus range (d).
Figure 3. Acoustic fields of the pseudolinear-speed waveguide with a pressure-free bottom calculated by analytical solution (a), KRAKEN (b), and SEM (c); TL curves at a depth of z = 40 m versus range (d).
Jmse 12 01499 g003
Figure 4. Sound speed of the waveguide with two-layer media and Munk waveguide (the blue line represents the sound speed profile in the media).
Figure 4. Sound speed of the waveguide with two-layer media and Munk waveguide (the blue line represents the sound speed profile in the media).
Jmse 12 01499 g004
Figure 5. Mode shapes for the waveguide with two-layer media, as calculated by KRAKENC (blue lines) and SEM (red dashed lines) at 33 Hz.
Figure 5. Mode shapes for the waveguide with two-layer media, as calculated by KRAKENC (blue lines) and SEM (red dashed lines) at 33 Hz.
Jmse 12 01499 g005
Figure 6. Acoustic fields of the waveguide with two-layer media and a pressure-free bottom calculated by KRAKEN (a), COUPLE (b), and SEM (c); TL curves at depths of z = 100 m versus range (d).
Figure 6. Acoustic fields of the waveguide with two-layer media and a pressure-free bottom calculated by KRAKEN (a), COUPLE (b), and SEM (c); TL curves at depths of z = 100 m versus range (d).
Jmse 12 01499 g006
Figure 7. Acoustic fields of the Munk waveguide with an acoustic half-space calculated by KRAKENC (a), COUPLE (b), and SEM (c); TL curves at a depth of z = 1000 m versus range (d).
Figure 7. Acoustic fields of the Munk waveguide with an acoustic half-space calculated by KRAKENC (a), COUPLE (b), and SEM (c); TL curves at a depth of z = 1000 m versus range (d).
Jmse 12 01499 g007
Figure 8. Horizontal wavenumbers for Munk waveguide, as calculated by KRAKENC and SEM.
Figure 8. Horizontal wavenumbers for Munk waveguide, as calculated by KRAKENC and SEM.
Jmse 12 01499 g008
Figure 9. Mode shapes for Munk waveguide, as calculated by KRAKENC (blue lines) and SEM (red dashed lines).
Figure 9. Mode shapes for Munk waveguide, as calculated by KRAKENC (blue lines) and SEM (red dashed lines).
Jmse 12 01499 g009
Table 1. The value of k r in KRAKEN and SEM.
Table 1. The value of k r in KRAKEN and SEM.
Numerical Models
KRAKEN SEM
Mode 10.2047404199838640.204740427046646
Mode 20.1971114277839660.197111422377371
Mode 30.1841926574707030.184192655388577
Mode 40.1643893569707870.164389353828404
Mode 50.1347005218267440.134700528071927
Mode 60.0853780359029770.085378039545094
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zhang, Y.; Tu, H.; Wang, Y.; Xu, G.; Gao, D. A Normal Mode Model Based on the Spectral Element Method for Simulating Horizontally Layered Acoustic Waveguides. J. Mar. Sci. Eng. 2024, 12, 1499. https://doi.org/10.3390/jmse12091499

AMA Style

Zhang Y, Tu H, Wang Y, Xu G, Gao D. A Normal Mode Model Based on the Spectral Element Method for Simulating Horizontally Layered Acoustic Waveguides. Journal of Marine Science and Engineering. 2024; 12(9):1499. https://doi.org/10.3390/jmse12091499

Chicago/Turabian Style

Zhang, Yinuo, Houwang Tu, Yongxian Wang, Guojun Xu, and Dongbao Gao. 2024. "A Normal Mode Model Based on the Spectral Element Method for Simulating Horizontally Layered Acoustic Waveguides" Journal of Marine Science and Engineering 12, no. 9: 1499. https://doi.org/10.3390/jmse12091499

APA Style

Zhang, Y., Tu, H., Wang, Y., Xu, G., & Gao, D. (2024). A Normal Mode Model Based on the Spectral Element Method for Simulating Horizontally Layered Acoustic Waveguides. Journal of Marine Science and Engineering, 12(9), 1499. https://doi.org/10.3390/jmse12091499

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop