Direct FVM Simulation for Sound Propagation in an Ideal Wedge

The sound propagation in a wedge-shaped waveguide with perfectly reflecting boundaries is one of the few range-dependent problems with an analytical solution. This provides a benchmark for the theoretical and computational studies on the simulation of ocean acoustic applications. We present a direct finite volume method (FVM) simulation for the ideal wedge problem, and both time and frequency domain results are analyzed. We also study the broadband problem with large-scale parallel simulations. The results presented in this paper validate the accuracy of the numerical techniques and show that the direct FVM simulation could be applied to large-scale complex acoustic applications with a high performance computing platform.


Introduction
The research on the sound propagation problem has a multidisciplinary and practical significance over oceanology, biology, shipbuilding, military affairs, and many other science and engineering subjects [1][2][3].Field survey and physical experiments are main approaches to the research of the underwater sound propagation.However, the high costs, low data rates, and difficulties in arranging experiments restrict the applications of experimental approach.With the rapid development of high performance computing technology, computational solution has become a much more important approach for the studying of underwater sound propagation problems.
The construction of a computational solution for this problem has been extensively investigated over the last decades.Underwater sound propagation can be mathematically described by the wave equation; thus, essentially different computational models use different approximations or discretization methods.In overall, there are seven types of computational methods to solve the sound propagation problem: the spectral or fast field program (FFP) [4], the normalmode solution (NM) [5][6][7], the parabolic equation solution (PE) [8][9][10], the ray modeling and ray solution [11][12][13], the wave-number integration method [14], and finite-difference method (FDM) [15,16] or finite element method (FEM) [17,18].In order to improve the computational efficiency, most of the methods (except FDM and FEM) reduce the wave equation into an approximate one.However, the efficiency of this type of approach is obtained by sacrificing generality through the applied assumptions and approximations.Reference [19] has had an in-depth discussion over the capabilities and limitations of those sound propagation models.
Running benchmark problems and comparing the numerical results are important and necessary steps for validating a computational solution.The ideal wedge problem describes the sound propagation in a wedge-shape waveguide, and it is one of the benchmarks proposed by Acoustical Society America (ASA) [20,21].It is one of the few range-dependent problems with an analytical solution [22,23].Hence, the ideal wedge problem has been extensively studied for validating the solutions calculated by various numerical approaches [18,24,25].Luo et al. [26] completed the original benchmark problem and provided an analytical solution for an ideal wedge problem with a rigid or a pressure-release bottom.
Although these fast computational methods were successfully used in explaining many observed ocean acoustic phenomena, there are still lots of scientific issues which cannot be accurately addressed by the approximated approach.For realistic problems, we usually have to solve a two-way wave equation with complex geometry.For this purpose, the numerical approach which directly discretizes the governing equations should be adopted.Direct FDM and FEM simulations [15,17] are used to model the acoustic problems.Nevertheless, since the direct discrete solutions must be able to represent the actual spatial and temporal variation of the acoustic field, all these methods are much more computationally intensive compared to the approximation methods.Therefore, the direct methods are rarely used to simulating the ocean acoustic propagation applications in the past.
Over the past few decades, high performance computing (HPC) techniques have made great achievements and the simulation time is possible to be significantly reduced through large-scale parallel computing, which means that the traditional direct simulation approach may become applicable for realistic ocean acoustic problems on a modern HPC platform.
In this paper, we design an acoustic numerical solver based on the finite volume method (FVM) and present the full numerical results of the direct FVM simulations for sound propagation in an ideal wedge.Both time domain and frequency domain results are analyzed; furthermore, we also model the broadband problem.The broadband probem is relatively difficult to simulate through the approximated approach; therefore, it is seldom studied over the last decades.The method proposed in this paper has its advantages for complex applications such as irregular boundary problems or broadband problems.
The remainder of the paper is organized as follows: the numerical techniques used in this paper are introduced in Section 2. Also a presentation of the analytical solution to the ideal wedge problem from the Acoustical Society of America (ASA) is included in Section 2. The 2D numerical results are presented and analyzed in Section 3. Accuracy analysis and large-scale parallel simulations are provided too.The conclusion follows in Section 4.

Model and Numerical Method
In this section, we give a synopsis of the ASA wedge benchmark and an introduction to the numerical techniques used in our simulations.

The Ideal Wedge Problem.
In the ideal wedge problem investigated by this paper, a pressure-release sea surface and a rigid or pressure-release bottom constitute the main boundaries of the 2D wedge.Both the surface and the bottom are perfectly reflecting boundaries and in the opposite the left boundary has to be set to an nonreflecting (absorbing) boundary.Although, in the real ocean environment, the wedge angle is very small, in this paper, we set it to 45 ∘ for simplicity.As shown in Figure 1, there is a line source situated in the wedge.
The 2D wedge problem is a vertical section of the 3D problem.The 3D wedge problem may have two types of wave sources.One is a spatial point locating in the field; the other takes a form of horizontal infinite straight line across the section parallel to the apex [26].The latter is an ideal, theoretical model not describing any entities directly.However, it plays an elemental and significant role like the 2D airfoil in flight vehicle aerodynamics and is widely used for experimental researches in computational acoustics.According to Luo et al. [5], even in 2-dimensional simulation, the problem with a point source generally has a cylindrical geometry, although the problem with a line source is usually described with cartesian coordinates.The solutions of these two problems take different forms.For the convenience of setting up the coordinates and meshing the computing domain, we choose the line source.
Luo et al. [26] proposed an analytical solution for the ideal wedge problem which can be used as a reference for accuracy validating.We summarized the solution in this section.
Under the 2D cylindrical coordinates, the ideal wedge problem could be expressed as where r = (,), ,  and  denote wave speed, mass density, and sound source.Applying the method of separation of variables, the wave equation could be reduced to the Helmholtz equation: The analytical solution to the Helmholtz equation with rigid and/or pressure-release bottom boundary conditions becomes Shock and Vibration 3 in which In order to obtain the exact value of the pressure level at each spatial point, this analytical solution also needs some numerical computations for the expansion coefficients of the Bessel and Hankel functions [23].

Numerical Technique.
The normal-mode methods and other frequency domain methods focus on the solution of the Helmholtz equation.These approaches solve the sound propagation process in the frequency domain.Nevertheless, through the direct FVM simulations, we can directly obtain the computational solution of the original wave equation.
To numerically solve the wave equation, we use an open source CFD toolbox released by the OpenCFD Ltd., named OpenFOAM.The equation is discretized through the finite volume method which locally satisfied the physical conservation laws through the integral over a control volume.For the spatial discretization terms, various predefined schemes are selectively applied and the temporal terms are discretized using a simple Euler scheme.Finally, the wave equation is reduced to a linear system; thus, using the iterative solvers in OpenFOAM, we can get the solutions for the equations at every time step.Solvers in the toolbox include the conjugate (PCG) and biconjugate gradient (PBiCG) methods.More details are presented in the OpenFOAM Manual.The boundary conditions for sea surface, bottom, and the left vertical numerical boundary are described numerically [27]: (i) The rigid bottom: (ii) The pressure-release bottom: (iii) The nonreflecting boundary: The numerical schemes used to descretize the wave equation and the boundary conditions are listed in Table 1.The solver and preconditioner used for solving  from the linear algebraic equations are listed in Table 2.

Results
In the following section, we first give a description to the platform and the parameters of the numerical experiments.Then, we give the numerical results for the periodic single frequency source problem with different bottom boundaries.
The problem with a broadband pulse source is also considered and tested.Finally, the large-scale parallel simulation is applied and analyzed.

Platform and Parameters
The Platform.In the experiments, we use an HPC cluster situated in the State Key Laboratory of High Performance Computing.This computing platform consists of hundreds of computing nodes and each node contains 12 Intel Xeon 2.1 GHz E5-2620 CPU cores and a total memory of 16 GB.
The Parameters.The parameters include the sound speed in the waveguide  = 1500 m/s, the density of the media  = 1.0 g/m 3 , the wedge angle of 45 ∘ , and the source location of 400 m in range and 200 m in depth.The periodic source problem uses a source with 25 Hz frequency.The broadband problem employs a square pulse lasting 0.1 s.

Accuracy
Analysis.The accuracy of the numerical techniques used in this paper is measured with the ASA ideal wedge problem.The specific features of this problem are shown in Figure 1.

Convergence Analysis.
Generally, the accuracy of the numerical simulation is dependent on its mesh density.A convergence analysis with different mesh desities is illustrated in Figure 2. In this section, we use three sets of mesh data, for a total of 127200, 508800, and 2035200 mesh cells, respectively.The original data are collected from the numerical pressure field  at 150-m depth as shown in Figure 2 where the  is the location of a field point, and the reference pressure takes the form of In the equation of  0 (),  0 is defined as  0 = 21500 with  denoting the source frequency [28].
The analytical solution to the Helmholtz equation is a frequency domain function.The direct FVM simulation outputs the time series of the pressure field  at every cell center of the mesh.In order to compare the results of the FVM simulation to the analytical solution, we take the Fourier transform of the FVM results from  = 2.0 s to  = 5.0 s, and we consider 2 receiver lines in depth of 30 m and 150 m.The FVM results are also shown in terms of TL defined as where the reference pressure  0 is set to the same value with the amplitude of the source wave.The locations of ridges and troughs of the FVM results are same with the analytical solution and remain invariant while the number of cells is increasing.This analysis shows that a mesh with 508800 cells is sufficient for the convergence.
The troughs of the curves in 2 indicate the destructive interference of the wave coming from the source and each reflected waves.Looking into the details, we may find that some numerical troughs (i.e., the one marked with an "X" in Figure 2(c)) go deeper than the analytical solution, while some others (i.e., the one marked with a "Δ" in Figure 2(c)) stop at a high position.And some troughs (i.e., the ones marked with "X"s in Figures 2(b) and 2(c)) go deeper when the field is divided into more cells.
Comparing those figures with Figure 3(f), it is clear that the troughs are horizontally located at the root of the pressure-range function.That is because (10) maps the () near 0 to minus infinity.This becomes the origin of the deep dips of the curves and of the differences of the numerical results and the analytical solution mentioned above.

Results for the Rigid Bottom.
Here, we consider the ideal wedge problem drawn in Figure 1 with a rigid bottom.Two snapshots of the complete spatial result of the direct FVM simulation for the rigid bottom problem are plotted in Figure 3.We also plot the the pressure field  at the time of 0.1 s, 1.0 s, and 2.0 s over the line of 150 m in depth.
The asymptotic stability of the FVM result is evident in Figure 3.In the initial 1 s time interval, the numerical solution presents the propagation procedure of the sound wave.The high similarity of the results in snapshot of  = 1.0 s and  = 2.0 s indicates the wave has reached its stable status in the waveguide.This fact rationalizes the usage of the time interval from  = 2.0 s to  = 5.0 s in the Fourier transform of the FVM results.
In Figure 4, the TL results in different depth are illustrated.The Fourier transform of multiple section-results could outline the complete features of the  field both spatially and temporally.Figures 4(a) and 4(b) show the comparison at 30-m depth and 150-m depth.For both the Fourier transform and the analytical solution, their TL forms are defined as that in ( 8) and (10).
Over both of the two receiver lines, the locations of ridges and troughs in range are the same between the FVM result and the analytical solution.Also, the error of the FVM simulation at most ridges has been controlled within 1 dB.

Results for the Pressure-Release Bottom.
In this experiment, we consider the ideal wedge problem with a pressurerelease bottom.Two snapshots of the spatial result of the direct FVM simulation to the problem are shown in Figure 5. Figure 5 also shows the spatial section of the pressure field  at the time of 0.1 s, 1.0 s, and 2.0 s along the line of 150 m in depth.
The results of the pressure-release bottom problem also show a behavior of temporal convergence.The wave propagation has reached the stable status after the time of 1.0 s.Observing the curves' right ends, we find a different reflection mode with that in the simulation for the rigid bottom problem.
In Figure 6, the TL results of the pressure-release bottom problem are illustrated.The results are selected at 30-m and 150-m depth.The Fourier transform of multiple sectionresults outline the complete features of the  field both spatially and temporally.Figures 6(a From the results shown in Figure 6, we can see that mostly the horizontal locations of the ridges and troughs in the FVM results meet a good agreement with those of the analytical solution.The errors at the ridges are controlled well.

Large-Scale Parallel Simulations.
In this section, the method of the direct FVM simulation is applied to the broadband problem which is seldom studied in the research on the approximated approach during the past few dacades.The ideal wedge waveguide is also used, and the results are analyzed after the experiments.between the sound speed  and the waveguide's span, the snapshots are selected from  = 0.03 s to  = 0.24 s.
The tone scale of the figures indicates the sound pressure .We can see that the sound wave has become flatter with its expansion.The width of the wave packet becomes larger over time.Intuitively speaking, the numerical simulation performs well on portraying the dissipation and dispersion of the sound wave in its propagation.

Conclusion
In this paper, we propose the platform and techniques of a direct FVM simulation for the sound propagation problem.The techniques of this simulation are applied to an ideal wedge problem characterized by a homogeneous water column, a pressure-release sea surface, a rigid or pressure-release bottom, and a periodic single frequency or a broadband pulse line source.The accuracy of the simulation is analyzed by comparing the numerical results with an analytical solution.A series of parallel computing experiments on solving the broadband problem is implemented.
We have studied the time and frequency domain results of the FVM simulation for the ideal wedge problem and have compared the results with the analytical solution proposed by Luo et al. [26].These comparisons reveal the accuracy of the method proposed in this paper.The numerical results obtain a good agreement with the analytical solution.The mesh convergence of the FVM simulation is analyzed.The convergence tests reveal that the mesh density does not have notable influence on the accuracy of the simulation.
Both the periodic single frequency source and the broadband pulse source are simulated and analyzed.We have investigated the time and frequency domain results.The platform and techniques of the simulation show a huge adaptability to these two different sources.By observing the physical phenomena, the asymptotic stability of the solutions to the periodic problem is considered.This fact helps us ensure the integration interval of the Fourier transform in the accuracy analysis.The numerical pulse generated by and propagated in the simulation for the broadband problem behaves well.
The large-scale parallel tests are implemented.The results on up to 384 processors show a good scalability.The techniques of the direct FVM simulation could be easily applied to a parallel environment with hundreds of processors, and this application could significantly reduce the CPU time in the simulation.
In summary, a direct FVM simulation for the ideal wedge problem with a homogeneous wedge column, a pressurerelease surface, a rigid or pressure-release bottom, and a periodic or single pulse line source is proposed.With the experiments and analysis, this simulation method shows its application prospects in the complicated simulations on ocean acoustics with higher performance computing platforms.

3 P r e s s u r e -r e l e a s e o r r i g i d b o t t o mFigure 1 :
Figure 1: The sketch of the ideal wedge environment used in the numerical experiments.The line source is located at 400 m in range and 200 m in depth.The source frequency is 25 Hz.

Figure 2 :
Figure 2: Comparison of the Fourier transform of the FVM result  at 150-m depth and the section of the analytical solution to the Helmholtz equation at the same depth: (a) The section line; (b) TL result in a mesh of 127200 cells; (c) TL result in a mesh of 508800 cells; (d) TL result in a mesh of 2035200 cells.The continuous line indicates the analytical solution and the descrete dots indicate the result of the FVM simulation.

Figure 3 :
Figure 3: The overview of the pressure field  of the rigid bottom problem at time 0.1 s and 1.0 s and the spatial section of the pressure field  at depth of 150 m: (a) and (d) show the snapshot at  = 0.1 s; (b) and (e) show the snapshot at  = 1.0 s; (c) and (f) show the snapshot at  = 2.0 s.
) and6(b)  show the comparison at 30-m depth and 150-m depth.

Figure 4 :Figure 5 :
Figure 4: Comparisons of the Fourier transform of the FVM simulation and the analytical solution to the rigid bottom problem: (a) TL result at depth of 30 m; (b) TL result at depth of 150 m.

3. 3 . 1 .Figure 6 :
Figure 6: Comparisons of the Fourier transform of the FVM simulation and the analytical solution to the pressure-release bottom problem: (a) TL result at depth of 30 m; (b) TL result at depth of 150 m.

Figure 7 :
Figure 7: The simulation for the single pulse source problem in the ideal wedge geometry.The bottom boundary is set rigid.The t value of each snapshot is listed in the subtitles of the figures.The small blue rings around the source in (g) and (h) are the secondary reflections on the edge of the source.

Table 1 :
Numerical schemes used for discretizing the wave equation and the boundary conditions.

Table 2 :
Solver and preconditioner used for solving  of the linear algebraic equations.