Two Chebyshev Spectral Methods for Solving Normal Modes in Atmospheric Acoustics
aa r X i v : . [ c s . C E ] J a n Two Chebyshev Spectral Methods for SolvingNormal Modes in Atmospheric Acoustics
Tu Houwang, a Wang Yongxian, b Xiao Wenbin, Lan Qiang, and Liu Wei College of Meteorology and Oceanography, National University of Defense Technology, Changsha,410073, China (Dated: 18 January 2021)
The normal mode model is important in computational atmospheric acoustics. It is oftenused to compute the atmospheric acoustic field under a harmonic point source. Its solutionconsists of a set of discrete modes radiating into the upper atmosphere, usually related to thecontinuous spectrum. In this article, we present two spectral methods, the Chebyshev–Tauand Chebyshev–Collocation methods, to solve for the atmospheric acoustic normal modes,and corresponding programs were developed. The two spectral methods successfully trans-form the problem of searching for the modal wavenumbers in the complex plane into a simpledense matrix eigenvalue problem by projecting the governing equation onto a set of orthog-onal bases, which can be easily solved through linear algebra methods. After obtaining theeigenvalues and eigenvectors, the horizontal wavenumbers and their corresponding modescan be obtained with simple processing. Numerical experiments were examined for bothdownwind and upwind conditions to verify the effectiveness of the methods. The runningtime data indicated that both spectral methods proposed in this article are faster than theLegendre–Galerkin spectral method proposed previously. © [https://doi.org(DOI number)][XYZ] Pages: 1–10 Keywords:
Chebyshev polynomial; normal modes; Taumethod; collocation method; computational atmosphericacoustics.
I. INTRODUCTION
The propagation of sound waves in the atmosphere isa basic subject of atmospheric acoustics . Sound waves inthe atmosphere undergo a series of complex processes, in-cluding ground reflection, atmospheric scattering, refrac-tion, and absorption . In fact, the propagation of soundwaves in the atmosphere satisfies the wave equation, butit is difficult to strictly solve the wave equation. Thus,scientists make approximations to the wave equation forspecific situations, thereby obtaining easy-to-solve equa-tions, which can be solved numerically to obtain a so-lution of the sound field. Numerical sound fields havethe advantages of intuitiveness and clarity, and they arewidely used in acoustic research. Based on this idea ofsolving the numerical sound field, computational atmo-spheric acoustics, a sub-discipline of atmospheric acous-tics, has been developed. Numerical models have manyforms. Different models are suitable for different environ-ments, and the results are not exactly the same. Main-stream numerical models include the parabolic equation(PE) model, the wavenumber integration method (the a [email protected] b [email protected] fast field program (FFP)) , and ray and Gauss beam approaches. The normal mode model is also a funda-mental method for solving for the acoustic field in theatmosphere with a finite ground impedance and horizon-tally stratified sound speed . A horizontally stratifiedatmosphere allows the wave equation to be solved by theseparation of variables method. After using Hankel in-tegral transforms, the sound field can be expressed interms of the sum of normal modes. When the groundimpedance is complex or there is sound attenuation inthe atmosphere, it is complicated to use the finite differ-ence method to solve for the atmospheric normal modes,and the result is not very accurate .In recent years, progress has been made on usingspectral methods to solve underwater acoustic problems,and small-scale research has begun to link the spectralmethods with the normal modes of underwater acous-tics. Dzieciuch developed MATLAB code for comput-ing normal modes based on Chebyshev approximations.Although he only realized the calculation of the simpleMunk waveguide, this was the first step in the appli-cation of the spectral methods to computational oceanacoustics. In 2016, Evans used the Legendre–Galerkinspectral method to develop a sound propagation calcu-lation program in a layered ocean environment. Sub-sequently, Tu et al. used the Chebyshev–Tau spec-tral method to develop a program for calculating soundpropagation in single-layer and layered ocean environ-ments. They subsequently solved for the normal modes inunderwater acoustics using the Chebyshev–Collocationmethod and proved that both of the spectral methods J. Acoust. Soc. Am. / 18 January 2021 1 ad high accuracy . They also applied the spectralmethods to the parabolic approximation of underwateracoustics . The results of these studies indicatedthat it is feasible to apply spectral methods for the cal-culation of underwater acoustics, and in many cases, ithas higher accuracy than the finite difference method.Monographs on spectral methods have also confirmedthis . Throughout the history of the developmentof atmospheric acoustics, many methods in underwateracoustics have been introduced . In computational at-mospheric acoustics, spectral methods are rarely used tocalculate the numerical sound field. In 2017, Evans successfully introduced the Legendre–Galerkin spectralmethod to construct atmospheric acoustic normal modes.He then further improved the method and proved theconvergence of the method .In this article, we propose two spectral methods forcalculating atmospheric acoustic normal modes. The re-sults are compared with Evans’s code , the correctnessof the two spectral methods proposed in this article wasverified, and computational speeds of the two spectralmethods were demonstrated. The text is organized asfollows. Section 2 describes normal modes in the atmo-sphere mathematically. Section 3 provides brief descrip-tions of the Chebyshev–Tau and Chebyshev–Collocationspectral methods and introduces the discretization of at-mospheric acoustic normal modes. In Section 4, two nu-merical experiments are shown to verify the correctnessof the methods proposed in this article. Section 5 an-alyzes the running speed of the spectral methods, andSection 6 concludes this article. II. ATMOSPHERIC NORMAL MODES
Acoustic theory reveals that the core of solving theacoustic field with a time-independent harmonic pointsource is the following wave equation : ρ ∇ · (cid:18) ρ ∇ p (cid:19) + k p = 0 . (1)In the above homogeneous Helmholtz equation, ρ is thedensity of the media, p is the sound pressure in thefrequency domain to be solved, and k represents thewavenumber, which is related to the frequency of thesource and the spatial position. ρ , p , and k are all func-tions of the spatial position, i.e., ρ ( x, y, z ), p ( x, y, z ), and k ( x, y, z ), respectively.We consider the medium of sound propagation to bethe atmosphere depicted in Figure 1. Taking the ∇ op-erator in Eq. (1) in cylindrical coordinates to obtain theacoustic governing equation in the cylindrical coordinatesystem ( r, z ), where r is the range, and z is the depth.Considering the case in Figure 1 where the density ρ ( z )and wavenumber k ( z ) are only related to depth (range-independent), Eq. (1) becomes:1 r ∂∂r (cid:18) r ∂p∂r (cid:19) + ρ ( z ) ∂∂z (cid:18) ρ ( z ) ∂p∂z (cid:19) + k ( z ) = 0 , (2) c(z) (z) (z)AtmosphereArtificial absorber0hHz Ground ImpedanceAir Impedance FIG. 1. Atmospheric sound propagation environment. where k ( z ) = ω/c ( z ), ω = 2 πf is the angular fre-quency of the sound source, f is the frequency of thesource, and c ( z ) is the sound speed profile. When con-sidering the attenuation of sound waves by the atmo-sphere, k ( z ) = [1 + iηβ ( z )] ω/c ( z ), where β ( z ) is the at-tenuating coefficient in units of dB per wavelength, and η = (40 π log e ) − . Through separation of variables, theacoustic pressure p ( r, z ) can be decomposed as follows: p ( r, z ) = ψ ( z ) R ( r ) , (3)where R ( r ) can be approximated by an analytical formof the function, and ψ ( z ) satisfies the following modalequation: ρ ( z ) dd z (cid:18) ρ ( z ) d ψ ( z )d z (cid:19) + k ( z ) ψ ( z ) = k r ψ ( z ) . (4)The modal equation is a Sturm–Liouville equation, andits characteristics are well known, that is, after addingappropriate boundary conditions, it has a series of modalsolutions ( k r , ψ ), where k r is a constant. When theconsidered medium has attenuation, k ( z ) is a complexfunction. The lower boundary of the atmosphere is theground, and sound waves on the ground usually need tomeet the following impedance boundary conditions:d ψ (0)d z + Gψ (0) = 0 , G = ik (0) /Z, (5)where Z is the normalized ground impedance. The up-per boundary of the atmosphere can be regarded as a freeboundary at infinity, or it can be called an acoustic half-space condition. To make the problem finite and solvablevia spectral methods, we add an artificial absorber layer[ h, H ] above the interest area [0 , h ], where the acousticparameters in [ h, H ] and [0 , h ] must be continuous. Theabsorber layer is usually set to be thick enough to at-tenuate the sound energy propagating upward, and noenergy is reflected back to the area of [0 , h ]. In this way,the following air impedance condition should be satisfiedat z = H :d ψ ( H )d z + αψ ( H ) = 0 , α = − ik ( H ) . (6)Solving the standard Sturm–Liouville problem will yieldmultiple sets of solutions ( k r m , ψ m ) , m = 1 , , · · · , where r m is called the horizontal wavenumber, and ψ m is calledthe eigenmode or mode. The modes of Eq. (4) are arbi-trary up to a nonzero scaling constant, so they should benormalized as follows: Z H [ ψ m ( z )] ρ ( z ) d z = 1 , m = 1 , , . . . . (7)Finally, the fundamental solution to the acoustic govern-ing equation (2) in the atmosphere can be approximatedas follows : p ( r, z ) ≈ √ π M X m =1 ψ m ( z s ) ψ m ( z ) exp( k r m r ) p k r m r , (8)where M is the number of modes used to synthesize thesound field.The core of solving for the normal modes of atmo-spheric acoustics is the solution of the differential equa-tions in Eq. (4)–(6). Solving for the normal modes ofthe atmospheric acoustics requires the discretization ofEq. (4)–(6). Traditionally, the domain of the problemsolved by the spectral method is usually in the interval[ − , x = 2 z/H − z ∈ [0 , H ] x ∈ [ − , z/ d x = H/
2, welet the operators L , P , Q have the following forms: L = (cid:20) H ρ ( x ) dd x (cid:18) ρ ( x ) dd x (cid:19) + k ( x ) (cid:21) , P = (cid:20) H dd x + G (cid:21) , Q = (cid:20) H dd x + α (cid:21) . Eq. (4)–(6) can be written in the following form: L ψ ( x ) = k r ψ ( x ) , x ∈ ( − , , P ψ ( −
1) = 0 , Q ψ (1) = 0 . (9)Next, we will develop two spectral methods to solve thissystem. III. DISCRETIZED ATMOSPHERIC NORMAL MODES BYTWO SPECTRAL METHODS
A spectral method is a kind of weighted residualmethod, and it can provide accurate solutions to differ-ential equations . In the spectral method, the un-known function to be solved ψ ( x ) is expanded by a setof linearly independent bases φ k ( x ). When the numberof bases tends to infinity, an accurate representation of ψ ( x ) can be obtained. However, in actual calculations,it is usually necessary to truncate to the first N -orderterms, thus obtaining an approximation of ψ ( x ), as fol-lows: ψ ( x ) = ∞ X k =0 ˆ ψ k φ k ( x ) ≈ ψ N ( x ) = N X k =0 ˆ ψ k φ k ( x ) , (10)where ˆ ψ k represents the expansion coefficients. Obtain-ing the value of ˆ ψ k is equivalent to obtaining the ap-proximate solution ψ N ( x ) of ψ ( x ). Inserting ψ N ( x ) from Eq. (10) into Eq. (9), Eq. (9) is no longer strictly true,and there is a residual Res ( x ), defined as follows: Res ( x ) = L ψ N ( x ) − k r ψ N ( x ) . (11)To make ψ N ( x ) as close to ψ ( x ) as possible, we needto minimize the residual through a certain principle .Setting the weighted integral of the residuals equal tozero is a widely used principle : Z − Res ( x ) w ( x )d x = 0 . (12)From Eq. (9), the residual Res ( x ) can be minimized onlyby adjusting the value of the expansion coefficients ˆ ψ k .The choice of the weight function w ( x ) is also crucial.In the two spectral methods developed in this article,the basis functions φ k ( x ) are both Chebyshev polyno-mials T k ( x ), and the difference is the selection of weightfunctions. The Chebyshev polynomial basis functions areprovided in the Appendix of this article. A. Discretized Atmospheric Normal Modes by Chebyshev–Tau Spectral Method
In the Chebyshev–Tau spectral method, in addi-tion to the basis functions being Chebyshev polynomials( φ k ( x ) = T k ( x )), the weight functions are also Chebyshevpolynomials ( w ( x ) = T k ( x )). Inserting Eq. (10) and (11)into Eq. (12), we obtain the new form of Eq. (12) for theChebyshev–Tau spectral method: Z − T j ( x ) √ − x L N X k =0 ˆ ψ k T k ( x ) − k r N X k =0 ˆ ψ k T k ( x ) ! d x = 0 j = 0 , , · · · , N − , .. (13)where √ − x is the orthogonal weighting factor of theChebyshev polynomial basis function space.This equation is also known as the weak form ofEq. (9). It will form ( N −
1) algebraic equations (ex-cluding the boundaries of x ), the two boundary condi-tions will produce two algebraic equations, and the un-knowns to be solved for are ˆ ψ to ˆ ψ N . The integral for-mulas listed in the above equations can be computed bythe Gauss–Lobatto quadrature to obtain accurate re-sults. To include the two end points of the domain x , theGauss–Lobatto nodes on x ∈ [ − ,
1] are taken : x j = − cos (cid:18) jπN (cid:19) , j = 0 , , · · · , N. (14) J. Acoust. Soc. Am. / 18 January 2021 3 here are two forms of Gauss–Lobatto quadrature : Z − f ( x ) √ − x d x ≈ N X j =0 f ( x j ) ω j , Z − f ( x )d x ≈ N X j =0 f ( x j ) ω j q − x j ,ω j = ( π N , j = 0 , N πN , otherwise , (15)where f ( x ) is the function to be integrated.We convert the original solution of the unknownfunction ψ ( x ) into solving for its expansion coefficientsˆ ψ k under the Chebyshev polynomial basis. The only dif-ficulty is the discretization of the operator L . The con-clusions used in the following text are directly given here.For detailed derivations, readers can refer to Refs. 17–20.The derivative dd x is included in the L operator, theexpanded coefficients ˆ ψ ′ k of ψ ′ ( x ) satisfy the followingrelationship with ˆ ψ k :ˆ ψ ′ k ≈ c k N X j = k +1 ,j + k =odd j ˆ ψ j , ˆ ψ ′ = ˆD ˆ ψ , (16)the second formula is the vector form of the first for-mula, where column vectors ˆ ψ = [ ˆ ψ , · · · , ˆ ψ N ] T and ˆ ψ ′ = [ ˆ ψ ′ , · · · , ˆ ψ ′ N ] T , respectively. ˆD is a square matrix oforder ( N + 1). To distinguish it from the differential ma-trix D in the Chebyshev–Collocation spectral method, ahat symbol is added to the relationship matrix.The known function k ( x ) is included in the L op-erator. Letting v ( x ) = k ( x ), there will be a productterm y ( x ) = v ( x ) ψ ( x ), and the expanded coefficients ˆ y k of y ( x ) satisfy the following relationship with ˆ ψ k :ˆ y k ≈ N X m + n = k ˆ ψ m ˆ v n + 12 N X | m − n | = k ˆ ψ m ˆ v n , ˆy ≈ ˆC v ˆ ψ . (17)Similarly, the second formula is the equivalent vectorform, ˆC v is also a square matrix of order ( N + 1), andthe subscript v indicates that the known function in theoperator is v ( x ).We show the discretization of the first term of theoperator L . We let g ( x ) = 1 ρ ( x ) , s ( x ) = g ( x ) ψ ′ ( x ) , f ( x ) = ρ ( x ) s ′ ( x ) . (18)In the Chebyshev–Tau spectral method, applyingEqs. (16) and (17) to Eq. (18), we can obtain ˆs = ˆC g ( ˆD ˆ ψ ) , ˆf = ˆC ρ ( ˆDˆs ) = ˆC ρ ( ˆD ( ˆC g ( ˆD ˆ ψ ))) . (19)The discrete forms of the operator L and Eq. (9) are asfollows: ˆL = (cid:18) H ˆC ρ ˆD ˆC g ˆD + ˆC v (cid:19) , ˆL ˆ ψ = k r ˆ ψ . (20) The boundary conditions produce algebraic equa-tions about the expansion coefficients in the Chebyshev–Tau spectral method as follows. In the Tau method, thefunction ψ ( x = ±
1) in the boundary conditions is also ex-panded by Eq. (10). The discretization of the boundaryoperators P and Q is similar to that of operator L , so thetwo boundary conditions generate two equations relatedto ˆ ψ k . To facilitate the description of the processing ofthe boundary conditions, the following intermediate rowvectors are defined as follows: ˆt = [ T ( − , T ( − , · · · , T N ( − , ˆp = 2 H ˆt ˆD + G ˆt , ˆt = [ T (1) , T (1) , · · · , T N (1)] , ˆq = 2 H ˆt ˆD + α ˆt . The matrix form of the discrete ground and air boundaryconditions Eq. (5) and (6) in the Chebyshev–Tau spectralmethod can be written as follows: ˆp ˆ ψ = 0 , ˆq ˆ ψ = 0 . (21)The algebraic equations formed by these two boundaryconditions and the ( N −
1) algebraic equations obtainedfrom the weak form are solved simultaneously, and wecan then solve for ˆ ψ k and obtain ψ ( x ).The row vectors ˆp and ˆq are used to replace the lasttwo rows of the ˆL matrix in Eq. (20), and the last twoelements of the column vector ˆ ψ on the right-hand sideof Eq. (20) are replaced with 0, so that the boundaryconditions are strictly met. We let the matrix composedof the first ( N −
1) rows and ( N −
1) columns of ˆL be ˆL .The matrix composed of the first ( N −
1) rows and thelast two columns of ˆL is ˆL . The row vectors composedof the first ( N −
1) elements of the row vectors ˆp , ˆq ,and ˆ ψ are ˆp , ˆq and ˆ ψ , respectively. The row vectorscomposed of the last two elements of the row vectors ˆp , ˆq , and ˆ ψ are ˆp , ˆq , and ˆ ψ , respectively. Thus, Eq. (20)can be changed to the following block form: ˆL ˆL ˆp ˆp ˆq ˆq ˆ ψ ˆ ψ N − ˆ ψ N = k r ˆ ψ . (22)According to the horizontal and vertical lines in the aboveformula, Eq. (22) can be abbreviated as follows: " ˆL ˆL ˆL ˆL ˆ ψ ˆ ψ = k r " ˆ ψ . (23)Eq. (23) can be solved as follows: ˆ ψ = − ˆL − ˆL ˆ ψ , ( ˆL − ˆL ˆL − ˆL ) ˆ ψ = k r ˆ ψ . (24)Therefore, a set of ( k r , ˆ ψ ) can be solved for by the( N − k r m , ˆ ψ m ), aneigensolution ( k r m , ψ m ) of Eq. (4) can be obtained byEq. (10). In this process, each eigenmode should be nor-malized by Eq. (7). Finally, the sound pressure field isobtained by applying Eq. (8) to the chosen modes. . Discretized Atmospheric Normal Modes by Chebyshev–Collocation Spectral Method The Collocation method uses the Dirac function δ ( x )as the weight function in Eq. (12). The characteristicsof the δ ( x ) function are well known. In the Collocationmethod, Eq. (12) becomes the following: Z − Res ( x ) δ ( x − x j )d x = Ri ( x j ) = L ψ ( x j ) − k r ψ ( x j ) ,j = 0 , , , ...N. (25)The above formula shows that in the Collocation method,the weighted residual principle becomes that the resid-uals are all 0 at the selected discrete points x j . Itsessence is to only make the original differential equation(9) strictly hold on this set of discrete points, so as tosolve for the function value ψ ( x j ) of the modal function ψ ( x ) on this set of discrete points as an approximation.In the Collocation method, there is no need to expandthe function to be sought as Eq. (10). This is why theCollocation method is considered to be a special spectralmethod, sometimes called the pseudospectral method .In the Chebyshev–Collocation method, we also take thediscrete points of the Chebyshev–Gauss–Lobatto nodesin Eq. (14). In this case, the only difficulty is the dis-cretization of operator L . The conclusions used in thefollowing text are directly given here as in the introduc-tion of Chebyshev–Tau spectral method. For a detailedderivation, readers can refer to Refs. 17–20.The derivative term ψ ′ ( x ) and ψ ( x ) have the follow-ing relationship: ψ ′ = D ψ , D = c j ( − j + l c k ( x j − x k ) , j = k − x k ( − x k ) , ≤ j = k ≤ N − N +16 , j = l = 0 − N +16 , j = l = Nc k = ( , k = 01 , k > , (26)where ψ ′ = [ ψ ′ ( x ) , ψ ′ ( x ) , ψ ′ ( x ) , · · · , ψ ′ ( x N )] T repre-sents the function value of the derivative term ψ ′ ( x ).Similarly, ψ = [ ψ ( x ) , ψ ( x ) , ψ ( x ) , · · · , ψ ( x N )] T . Ma-trix D is also called the Chebyshev–Collocation differen-tial matrix.The product y ( x ) = v ( x ) ψ ( x ) can be processed asfollows: y = C v ψ , (27)where y = [ y ( x ) , y ( x ) , y ( x ) , · · · , y ( x N )] T , C v is a ( N +1) × ( N + 1) diagonal matrix, and ( C v ) ii = v ( x i ) , i =0 , , ..., N .For the Collocation method, the boundary condi-tions are only related to the endpoints of the domain x , so the discrete points on the boundaries ( x and x N )only need to satisfy the boundary conditions, not the differential equation. The discretized forms of the oper-ators P and Q are similar to that of operator L . Sim-ilar to the Chebyshev–Tau spectral method, the oper-ator L also needs to be discretized in the Chebyshev–Collocation method. With reference to Eq. (18) and (19),in the Chebyshev–Collocation method, the L operatorand Eq. (9) have the following forms: L = (cid:18) H C ρ DC g D + C v (cid:19) , L ψ = k r ψ . (28)To facilitate the description of the processing of theboundary conditions, the first and last rows of D aredefined as row vectors d and d , respectively. The fol-lowing intermediate ( N + 1)-dimensional row vectors aredefined as follows: t = [1 , , · · · , , p = 2 H d + G t , t = [0 , , · · · , , q = 2 H d + α t . The matrix form of the discrete ground and airimpedance conditions Eq. (5) and (6) in the Chebyshev–Collocation method can be written as follows: p ψ = 0 , q ψ = 0 . (29)In the Collocation method, the row vectors p and q are used to replace the first row and the last row of the L matrix in Eq. (28), so that the boundary conditions aresatisfied. We let the block matrix formed by the secondrow to the N -th row of the matrix L be L , and thecolumn formed by the second to the N -th elements of ψ be ψ . Eq. (28) can then be written as follows: pL q ψ ψ ψ N = k r ψ . (30)We only need to perform a simple row transformationand column transformation on Eq. (30) to transformit to a form similar to Eq. (23), and then we use thesame method used for Eq. (24) to find the eigenval-ues/eigenvectors. A set of ( k r , ψ ) can be solved by the( N − IV. NUMERICAL EXPERIMENT AND ANALYSIS
To verify the correctness of the above presentedspectral methods in solving the normal modes ofatmospheric acoustics, the authors developed thecorresponding programs based on the above deriva-tion. The programs based on Chebyshev–Tau spectralmethod and Chebyshev–Collocation method arecalled ‘AtmosCTSM’ and ‘AtmosCCSM,’ respec-tively. The code was written in FORTRAN/MATLAB
J. Acoust. Soc. Am. / 18 January 2021 5 nd is available at the author’s GitHub homepage( https://github.com/tuhouwang/Atmospheric-normal-modes ).For comparison, we considered the program ‘aaLG’ basedon the Legendre–Galerkin spectral method, which wasdeveloped by Evans in FORTRAN and verified bycomparison with PE and FFP . The two examplesshown by Evans can be used as benchmark examples.The source frequency of both cases was 100 Hz at aheight of 5 m above the ground. The normalized groundimpedance Z (related to the constant G in Eq. (6)) isthe same as the value used by Gilbert , and the value is Z = 12 .
97 + 12 . i . In the following two experiments,the order of the spectral truncation N in the threespectral methods was taken as 1500. Using the TL toexpress the acoustic field , the relationship betweenit and the sound pressure is TL = −
20 log ( | p | / | p | ),where p is the sound pressure at a distance of 1 m fromthe sound source. A. Downwind Case
The first numerical experiment was a downwindcase. The piecewise linear acoustic parameter profileused in this numerical experiment, which was presentedby Evans , is shown in Table I. In contrast to the caseconsidered by Evans , the change of the atmosphericdensity with height was considered in this work. There-fore, a column of density data is included in the table.In fact, when the density is taken as a constant, ρ ( z )and 1 /ρ ( z ) in Eq. (5) will be eliminated, which meansthat the uniform density has no effect on the propaga-tion of normal modes. The table clearly reveals that thethickness of the atmosphere is 700 m, and the artificialabsorber is located between 700 and 2000 m. TABLE I. Piecewise linear acoustic parameter profile used inexperiment 1, cited from Evans . Height (m) Sound speed (m/s) Attenuation (dB/wavelength) Density (kg/m3)2000 344.0 2.50 constant1500 344.0 0.10900 344.0 0.01700 344.0 0.00500 341.5 0.00100 349.0 0.000 345.0 0.00
Figure 2 shows the horizontal wavenumbers k r cal-culated by the Legendre–Galerkin spectral method andthe two spectral methods developed in this article on thecomplex plane. The consistency of the eigenvalue dis-tribution in the figure illustrates the correctness of thehorizontal wavenumbers calculated by the three methods.Figure 3 shows the first four normal modes of experiment1. It reveals that the modes obtained by the two spectralmethods proposed in this article were highly consistentwith those obtained by the Legendre–Galerkin spectralmethod. Figure 4 presents an overview of the acous-tic fields obtained by the three methods. We used thefirst 552 modes with phase velocities between 341.7 and 391.2 m/s to synthesize the sound fields. The horizon-tal wavenumbers of these modes are shown in Figure 2.The acoustic fields calculated by the three methods werevery similar. Figure 5 shows the TL curves versus therange for a receiver at a height of 1 m over the rangeinterval 0–5 km. The results of the two spectral methodspresented in this article were very similar to those of theLegendre–Galerkin spectral method, and there may havebeen small differences only in the acoustic shadow areas. Real Wave Number (1/m) I m ag i na r y W a v e N u m be r ( / m ) Legendre-GalerkinChebyshev-TauChebyshev-Collocation
FIG. 2. Horizontal wavenumbers calculated by the Legendre–Galerkin spectral method and the two spectral methods pro-posed in this article for experiment 1. -0.2 0 0.202004006008001000 H e i gh t ( m ) Legendre-Galerkin Chebyshev-Tau Chebyshev-Collocation -0.2 0 0.202004006008001000 -0.2 0 0.202004006008001000 -0.2 0 0.202004006008001000 FIG. 3. First four normal modes calculated by the three spec-tral methods in experiment 1.
B. Upwind Case
The second numerical experiment is an upwind case.The piecewise linear acoustic parameter profile used inthis numerical experiment was presented by Evans , andit is shown in Table II. The table clearly illustrates thatthe thickness of the atmosphere was 900 m, and the ar-tificial absorber was located between 900 and 2000 m.Figure 6 shows four modes computed by theLegendre–Galerkin spectral method, Chebyshev–Tauspectral method, and Chebyshev–Collocation method. a)(b)(c) FIG. 4. Atmospheric acoustic fields obtained by (a)Legendre–Galerkin spectral method , (b) the Chebyshev–Tau spectral method, and (c) the Chebyshev–Collocationmethod of experiment 1. Range (km) T L ( d B ) Legendre-GalerkinChebyshev-TauChebyshev-Collocation
FIG. 5. TL versus range for a receiver at a height of 1 m overthe range interval 0–5 km from Legendre–Galerkin spectralmethod , Chebyshev–Tau spectral method, and Chebyshev–Collocation method.TABLE II. Piecewise linear acoustic parameter profile usedin experiment 2, cited from Evans . Height (m) Sound speed (m/s) Attenuation (dB/wavelength) Density (kg/m3)2000 346.0 1.00 constant1500 346.0 0.101200 346.0 0.01900 346.0 0.00500 348.0 0.00350 344.0 0.00100 340.0 0.000 344.0 0.00
The modes obtained by the three methods are drawnin the same figure. The three lines almost completelyoverlap in the subfigures, and the differences betweenthem are insignificant. Figure 7 presents the acousticfields calculated by the three spectral methods, where553 modes with phase velocities less than 393.2 m/s wereused to synthesize sound field. In the atmosphere layer,the sound fields calculated by the three methods werehighly consistent. Figure 8 shows the TL curves versusthe range for a receiver at a height of 1 m over the rangeinterval 0–10 km from the Legendre–Galerkin spectralmethod, the figure shows that the results of several meth-ods were very similar. The differences between the threemethods are indistinguishable at this plotting accuracy,and there may have been small differences only in theacoustic shadow area.From the numerical results displayed above, we cansee that three methods with different theoretical founda-tions all yielded very similar acoustic fields and normalmodes, regardless of whether the sound speed profile wasdownwind or upwind. The consistency of these threemethods proved that the two spectral methods proposedin this article are feasible for solving the atmospheric nor-mal modes.
J. Acoust. Soc. Am. / 18 January 2021 7 H e i gh t ( m ) Legendre-Galerkin Chebyshev-Tau Chebyshev-Collocation -0.2 0 0.2050100150200250 -0.2 0 0.2050100150200250 -0.2 0 0.2050100150200250 FIG. 6. (a)–(d) First, second, fourth, and sixth modes ob-tained by Legendre–Galerkin, Chebyshev–Tau spectral andChebyshev–Collocation methods for experiment 2.
V. DISCUSSION OF COMPUTATIONAL SPEED
To further compare the characteristics of the twospectral methods proposed in this article, we divided eachmethod into four steps, and we discuss the running timeand complexity of each part separately. The four stepsare as follows: discretizing the equation, solving eigen-value problems, obtaining normal modes, and synthesiz-ing the sound field. Table III lists the time consumptionof each step of the programs in the two experiments. Thetime listed in the table is the average of ten tests. In thetests, the three programs were run on a Dell XPS 8930desktop computer equipped with an Intel i7-8700K CPU.The FORTRAN compiler used in the test was gfortran7.5.0.In terms of speed, the AtmosCCSM was slightlyfaster than the AtmosCTSM. This is because the Taumethod requires forward and backward Chebyshev trans-formations, unlike the Collocation method. The aaLGprogram was much slower than the two newly devel-oped programs. Solving the eigenvalue problem was themost time-consuming step for the three programs. More-over, the aaLG program spent much more time thanthe other two programs on solving the eigenvalue prob-lem. However, the aaLG uses a subroutine developed byEvans to solve the matrix eigenvalue problem, whileboth the AtmosCTSM and AtmosCCSM solve the eigen-value problem by calling the Lapack numerical library. Infact, matrix eigenvalue problems have the same compu-tational complexity O ( N ). It is apparent from this tablethat the subroutine written by Evans is much slower thanthe Lapack numerical library, which is the main reasonthat the aaLG consumed much more time than the othertwo programs.We modified aaLG to also call the Lapack numericallibrary when solving the eigenvalue problem. The modi-fied program was named ‘aaLG-M’. The fourth column ofTable III lists the running time of aaLG-M. The aaLG-Mtook roughly the same time to solve the matrix eigenvalueproblem as the other two programs. However, the aaLG-M was still slower than the other two programs. Themost significant difference of the running time between (a)(b)(c) FIG. 7. Atmospheric acoustic fields obtained by (a)Legendre–Galerkin spectral method , (b) the Chebyshev–Tau spectral method, and (c) the Chebyshev–Collocationmethod for experiment 2. Range (km) T L ( d B ) Legendre-GalerkinChebyshv-TauChebyshev-Collocation
FIG. 8. TL versus range for a receiver at a height of 1m over the range interval 0–10 km from the (a) Legendre–Galerkin spectral method and the (b) Chebyshev–Tau spec-tral method and Chebyshev–Collocation method.TABLE III. Time consumption of each step of each programin two experiments in units of seconds. Experiment Part of program aaLG aaLG-M AtmosCTSM AtmosCCSMdownwind 1 105.344 104.691 0.522 0.4682 2017.324 34.289 34.091 34.3313 35.587 35.292 0.867 0.2374 10.021 8.714 0.518 0.421Total 2138.276 182.986 35.998 35.184upwind 1 125.429 123.892 0.482 0.3612 2039.324 36.119 34.886 34.0173 36.501 38.181 0.911 0.3344 11.669 10.648 0.806 0.616Total 2212.923 208.840 37.085 35.328 the three programs was in the first step (discretizing theequation). In the first step, each element of the ma-trix finally obtained by the Legendre–Galerkin spectralmethod must be numerically integrated for every piece-wise linear acoustic profile, which means the number ofcalculations is very large. In contrast, the two methodsproposed in this article only need to perform simple inter-polation of the acoustic profile and matrix multiplicationto obtain the discrete equations. In the mode obtain-ing and normalization steps, the two spectral methodsdevised in this article still required less time than theLegendre–Galerkin method. It is worth mentioning thatthe AtmosCTSM.m and AtmosCCSM.m programs (de-veloped in MATLAB, which is better at matrix opera-tions) could obtain the results of the above experimentsin less than 4 seconds (run on the same platform in MAT-LAB 2019a), which is an attractive result.
VI. CONCLUSION
In this article, we propose two spectral methods forsolving for atmospheric acoustic normal modes. An arti- ficial absorption layer was added above the atmosphereof interest to reduce the impact of the truncated half-space on the area of interest. Next, we designed twoexamples, performed a detailed analysis of the results ofeach example, and finally verified the correctness and re-liability of the proposed methods. Tests on the runningtime of the programs developed based on three methodsshowed that, in terms of the running time, the meth-ods proposed in this article had better speeds than theLegendre–Galerkin spectral method.
APPENDIX:
The Chebyshev polynomials are defined in the inter-val x ∈ [ − ,
1] and have the following definition: T ( x ) = 1; T ( x ) = x ; T ( x ) = 2 x − ,T k +1 ( x ) = 2 xT k ( x ) − T k − ( x ) , k = 2 , , , · · · . (A.1)The orthogonality of these polynomials is defined as fol-lows: Z − T j ( x ) T k ( x ) √ − x d x = , j = kπ, j = k = 0 π/ , j = k ≥ , (A.2)where √ − x is the orthogonal weighting factor. Theexpansion coefficients of a known function ψ ( x ) on theChebyshev polynomial basis can be obtained by the fol-lowing formula:ˆ ψ k = 2 πc k Z − ψ ( x ) T k ( x ) √ − x d x, c k = ( , k = 01 , k > . (A.3)Eq. (A.3) is called the forward Chebyshev transform, andit can be quickly calculated using the fast Fourier trans-form technique introduced by Canuto . Acknowledgments
The authors are very grateful to Richard B. Evansfor providing the aaLG program in Ref. 22.This study was supported by the National Key Re-search and Development Program of China [grant num-ber 2016YFC1401800] and the National Natural Sci-ence Foundation of China [grant numbers 61972406,51709267]. E. M. Salomons,
Computational atmospheric acoustics (SpringerScience Business Media, 2001). X. Yang,
Computational atmospheric acoustics (Science press,Beijing, 2015). K. Gilbert and M. White, “Application of the parabolic equationto sound propagation in a refracting atmosphere,” Journal of theAcoustical Society of America , 630–637 (1989). J. Acoust. Soc. Am. / 18 January 2021 9
K. E. Gilbert and X. Di, “A fast green’s function method forone-way sound propagation in the atmosphere,” The Journal ofthe Acoustical Society of America , 2343–2352 (1993). F. B. Jensen, W. A. Kuperman, M. B. Porter, and H. Schmidt,
Computational Ocean Acoustics (Springer-Verlag, New York,2011). R. Gilbert, S. W. Lee, E. Kuester, D. C. Chang, W. F. Richards,R. Gilbert, and N. Bong, “A fast-field program for sound prop-agation in a layered atmosphere above an impedance ground.,”Journal of the Acoustical Society of America , 345 (1985). M. B. Porter, “Scooter: A finite element ffp code” (2010), https://oalib-acoustics.org/AcousticsToolbox/index_at.html . A. D. Pierce,
Acoustics: An Introduction to its Physical Princi-ples and Applications (American Institute of Physics, New York,1991). M. Dzieciuch, “A matlab code for computing nor-mal modes based on chebyshev approximations” (1993), https://oalib-acoustics.org/Modes/aw . R. B. Evans, “A Legendre-Galerkin technique for differ-ential eigenvalue problems with complex and discontinu-ous coefficients, arising in underwater acoustics” (2016), https://oalib-acoustics.org/Modes/rimLG/ . H. Tu, Y. Wang, W. Liu, X. Ma, W. Xiao, and Q. Lan, “Achebyshev spectral method for normal mode and parabolic equa-tion models in underwater acoustics,” Mathematical Problems inEngineering , 1–12 (2020) doi: . H. Tu, Y. Wang, Q. Lan, W. Liu, W. Xiao, and S. Ma,“A chebyshev-tau spectral method for normal modes of un-derwater sound propagation with a layered marine environ-ment,” Journal of Sound and Vibration , 1–12 (2020) doi: . H. Tu, “A chebyshev-tau spectral method for normalmodes of underwater sound propagation with a lay-ered marine environment in matlab and fortran” (2021), https://oalib-acoustics.org/Modes/NM-CT . H. Tu, Y. Wang, Q. Lan, W. Liu, W. Xiao, and S. Ma, “Applyinga chebyshev collocation method based on domain decompositionfor calculating underwater sound propagation in a horizontallystratified environment,” (2020), Journal of Sound and Vibration(under review). Y. Wang, H. Tu, W. Liu, W. Xiao, and Q. Lan, “Application ofchebyshev collocation method to solve parabolic equation modelof underwater acoustic propagation,” (2020), Acoustics Australia(major revision). H. Tu, Y. Wang, X. Ma, and X. Zhu, “Applying chebyshev-tauspectral method to solve the parabolic equation model of wide-angle rational approximation in ocean acoustics,” (2020), Journalof Theoretical and Computational Acoustics (under review). D. Gottlieb and S. A. Orszag,
Numerical Analysis of SpectralMethods, Theory and Applications (Society for Industrial andApplied Mathematics, Philadelphia, 1977). J. P. Boyd,
Chebyshev and Fourier Spectral Methods (SecondEdition, Dover, New York, 2001). C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang,
Spec-tral Methods Fundamentals in Single Domains (Spring-Verlag,Berlin, 2006). S. Jie, T. Tao, and W. Lilian,
Spectral Methods Algorithms,Analysis and Applications (Springer-Verlag, Berlin, Heidelberg,2011). R. B. Evans, X. Di, and K. E. Gilbert, eds.,
A Legendre-Galerkintechnique for finding atmospheric acoustic normal modes , Vol. 30(Acoustical Society of America, Boston, Massachusetts, USA)(2017). R. B. Evans, X. Di, and K. E. Gilbert, “A legendre-galerkinspectral method for constructing atmospheric acoustic normalmodes,” The Journal of the Acoustical Society of America ,3595–3601 (2018). R. B. Evans, “The convergence of the legendre–galerkin spectralmethod for constructing atmospheric acoustic normal modes,” Journal of Theoretical and Computational Acoustics ,2050002 (2020) doi: . P. J. Anderson and G. Loizou,
A Jacobi type method for complexsymmetric matrices , Vol. 25 (Numerische Mathematik, 1976),pp. 347–363., Vol. 25 (Numerische Mathematik, 1976),pp. 347–363.