均勻海洋波導(dǎo)本征值問題的多模態(tài)分析方法

魏 旺，王 斌，范 軍

(上海交通大學(xué)高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，海洋工程國家重點(diǎn)實(shí)驗(yàn)室，上海 200240)

0 引 言

簡正波理論[1]是用來分析波導(dǎo)中聲傳播準(zhǔn)確解的常用方法，在單層等聲速波導(dǎo)中通常能夠得到本征值和本征函數(shù)的解析式從而簡化計(jì)算。但在一些復(fù)雜的情況下(例如具有聲速剖面的波導(dǎo)、雙層波導(dǎo))，由復(fù)雜邊界帶來的超越方程的搜根問題，使得本征值、本征函數(shù)的計(jì)算量很大。目前已有的常用的方法有圖像法、牛頓迭代法等，但這些方法的計(jì)算速度慢，精度低，且會(huì)受到設(shè)定初值的影響，從而產(chǎn)生丟根和漏根的情況。為了提高本征值的求解效率，譜方法[2]被應(yīng)用于求解波導(dǎo)中的超越方程，避免了傳統(tǒng)搜根法帶來的漏根，且具有較快的收斂速度，其中較為常用的是切比雪夫(Chebyshev)多項(xiàng)式進(jìn)行迭代計(jì)算，但該方法在某些參數(shù)選取不恰當(dāng)時(shí)帶來的收斂性并不理想[3]。

PAGNEUX等[4]、BI等[5]在計(jì)算聲管中的聲傳播時(shí)提出了多模態(tài)方法(Multimodal Approach,MA)，用正弦或者余弦函數(shù)作為譜方法中的正交基，求解出聲管中的本征函數(shù)，比搜根法具有更高的效率。雖然與切比雪夫多項(xiàng)式相比，多模態(tài)方法需要更高的截?cái)嚯A數(shù)，但其結(jié)果總是收斂的，在不同邊界條件下的均勻、非均勻波導(dǎo)中都具有一定的應(yīng)用前景。

本文將該方法推廣至海洋波導(dǎo)中，將超越方程的搜根問題轉(zhuǎn)化為正交展開系數(shù)矩陣的特征值分解。計(jì)算了單層及雙層波導(dǎo)的本征值與本征函數(shù)，通過與解析解以及簡正波程序KRAKEN的結(jié)果對(duì)比，驗(yàn)證了該方法的有效性，為進(jìn)一步研究均勻、非均勻波導(dǎo)中的聲傳播提供了新的思路。

1 均勻海洋波導(dǎo)中的多模態(tài)方法

如圖1所示的軸對(duì)稱均勻波導(dǎo)中，水平和深度方向分別是r軸和z軸，假設(shè)海面位于z=0處，海底位于z=h處。海底以下直到H是沉積層。波導(dǎo)中海水聲速cw(z)在豎直方向上連續(xù)變化，密度為常數(shù)ρw，沉積層的聲速和密度分別為常數(shù)cb和ρb。

圖1 均勻波導(dǎo)模型Fig.1 Model of the homogeneous waveguide

波動(dòng)方程的柱坐標(biāo)表達(dá)式為[1]

其中：k(z)=2πf/c(z)，c(z)在海水和海底中分別記為cw(z)和cb；f為聲源頻率；p代表聲壓，其滿足海面自由、海底處聲壓和法向振速的連續(xù)條件：

假設(shè)y=H處為剛性邊界，即

通常情況下，根據(jù)式(2)的前兩個(gè)條件，可以設(shè)本征函數(shù)具有式(4)的形式[6]

因此，本征方程的求解轉(zhuǎn)化為對(duì)超越方程(5)的搜根問題。一般的解法有圖像法、牛頓迭代法等，但這些方法速度慢、精度不高且依賴于初值。因此這里引入多模態(tài)展開方法[7]，令：

式中，Pm(r)為不同階的展開系數(shù)。正交函數(shù)ψm滿足

其中：

設(shè)mψ具有如式(4)所示的形式：

其中：

使用正弦函數(shù)的優(yōu)勢在于，其任意階導(dǎo)數(shù)都具有解析表達(dá)式，避免了數(shù)值發(fā)散(例如使用切比雪夫多項(xiàng)式[2])。計(jì)算時(shí)，通常取m=0,1,??,M-1的前M階進(jìn)行截?cái)?，?/p>

其中，Ψ(z)是由正交函數(shù)ψm(z)組成的M元列向量，上標(biāo)T表示轉(zhuǎn)置，P(r)表示不同階的展開系數(shù)Pm(r)組成的M元列向量。

其中，P′′(r)、P'(r)分別表示由組成的列向量，A是M×M維矩陣：

矩陣C、F、E中的元素分別為

矩陣A經(jīng)特征值分解得到

其中，Λ代表矩陣A的特征值dm組成的對(duì)角陣。X代表特征向量組成的列矩陣。

根據(jù)式(12)與式(17)可以得到

式(18)為零階貝塞爾方程，其具有通解形式

其中，+、-分別表示向外、向內(nèi)傳播的分量。D1和D2分別滿足：

對(duì)于二維波導(dǎo)情況，距離和深度方向分別為x和y軸，波動(dòng)方程滿足：

其通解為

其中：

關(guān)于本征值的超越方程搜根問題被轉(zhuǎn)換為矩陣的特征值問題，同時(shí)可以得到本征函數(shù)，因此尤其適合于本征函數(shù)無解析表達(dá)式的情況。

2 數(shù)值驗(yàn)證與分析

本節(jié)分別給出多模態(tài)方法在單層等聲速波導(dǎo)、具有聲速剖面的波導(dǎo)以及雙層波導(dǎo)的本征值與本征函數(shù)計(jì)算結(jié)果，并與一些標(biāo)準(zhǔn)結(jié)果進(jìn)行對(duì)比驗(yàn)證，其中通過本征函數(shù)(振型)的分布可以看出各個(gè)模態(tài)的傳播規(guī)律。

2.1 單層波導(dǎo)的本征值與本征函數(shù)

單層等聲速波導(dǎo)中，設(shè)水中密度ρw=1 000 kg·m-3，海底為水平邊界，水深H=500 m，分析頻率為25 Hz，等聲速情況下設(shè)其聲速cw=1 500 m·s-1為常數(shù)。由于是單層波導(dǎo)，令雙層模型中的ρb=ρw，ρb=ρw，即h=H，具體推導(dǎo)過程不再重復(fù)。

等聲速情況下，波導(dǎo)中的聲場分布情況可以寫成解析解的形式[8]，其中本征值可以表示為

表1給出了解析解與多模態(tài)展開方法的計(jì)算結(jié)果，由于單層波導(dǎo)中的本征方程及其導(dǎo)數(shù)有較好的連續(xù)性，因此用較低的模態(tài)階數(shù)即可得到較為準(zhǔn)確的計(jì)算結(jié)果，這里將截?cái)嚯A數(shù)選取為16。

表1 等聲速波導(dǎo)本征值Table 1 Eigenvalues of isovelocity waveguide

本征函數(shù)如圖2所示，可以看到多模態(tài)方法的本征值、本征函數(shù)計(jì)算結(jié)果與解析解具有很好的一致性。

圖2 等聲速波導(dǎo)本征函數(shù)Fig.2 Eigenfunctions in isovelocity waveguide

當(dāng)波導(dǎo)中隨深度方向具有聲速剖面時(shí)，假設(shè)該聲速剖面如圖3所示，聲速剖面的表達(dá)式為[9]

圖3 聲速剖面示意圖Fig.3 A typical sound velocity profile

在具有聲速剖面的情況下，本征值和本征函數(shù)沒有解析解，因此可以通過數(shù)值方法進(jìn)行計(jì)算，這里將KRAKEN計(jì)算程序的數(shù)值結(jié)果作為標(biāo)準(zhǔn)解與多模態(tài)展開方法進(jìn)行比較，波導(dǎo)中除聲速外其他參數(shù)與等聲速情況一致，其本征值計(jì)算結(jié)果如表2所示。在同樣的計(jì)算機(jī)運(yùn)行環(huán)境中，KRAKEN方法的計(jì)算時(shí)間為0.342 s，多模態(tài)展開算法的計(jì)算時(shí)間為0.017 s。另外，通過觀察本征值的數(shù)量級(jí)可以知道，相鄰模態(tài)的本征值之差很小，這就對(duì)搜根法提出了更為嚴(yán)格的要求。若搜根的步長過大，很容易產(chǎn)生漏根，若步長過小，則會(huì)產(chǎn)生較多重根，從而導(dǎo)致計(jì)算效率的降低。本征函數(shù)結(jié)果如圖4所示，可以看到兩種方法的計(jì)算結(jié)果相吻合。另外，通過對(duì)比圖2與圖4可得，等聲速情況下本征函數(shù)在深度方向上均勻分布，而在具有聲速剖面時(shí)，低階模態(tài)主要集中在低聲速區(qū)域傳播。

表2 變聲速波導(dǎo)本征值Table 2 Eigenvalues in the waveguide with sound velocity profile

圖4 變聲速波導(dǎo)本征函數(shù)Fig.4 Eigenfunctions in the waveguide with sound velocity profile

2.2 雙層波導(dǎo)的本征值與本征函數(shù)

下面討論雙層波導(dǎo)的本征值與本征函數(shù)。分層波導(dǎo)中的本征值與本征函數(shù)依然沒有解析解的形式。設(shè)總深度H=2 000 m，水深h=100 m，水中聲速為cw=1 500 m·s-1，密度為ρw=1 000 kg·m-3，海底沉積層中的聲速為cb=1 700 m·s-1，密度為ρb=1 000 kg·m-3。由于本征函數(shù)的導(dǎo)數(shù)在海底界面處的不連續(xù)性，通常需要較高的模態(tài)截?cái)嚯A數(shù)進(jìn)行計(jì)算以得到可靠結(jié)果。這里將模態(tài)的截?cái)嚯A數(shù)取為200，表3給出了兩種方法得到的本征值，其中KRAKEN計(jì)算時(shí)間為0.371 s，多模態(tài)展開算法的計(jì)算時(shí)間為0.061 s。本征函數(shù)結(jié)果如圖5所示，可以看到，不同的深度上各階模態(tài)的主要傳播區(qū)域有所不同，只有前2階模態(tài)主要在水中傳播。

表3 雙層波導(dǎo)本征值Table 3 Eigenvalues in double-layered waveguide

圖5 雙層波導(dǎo)本征函數(shù)Fig.5 Eigenfunctions in double-layered waveguide

3 結(jié) 論

本文提出了海洋波導(dǎo)中本征值、本征函數(shù)求解的多模態(tài)展開方法，將本征函數(shù)用正弦函數(shù)作為正交基進(jìn)行展開，將關(guān)于超越方程的搜根問題轉(zhuǎn)化為展開系數(shù)矩陣的本征值求解問題，從而避免了傳統(tǒng)方法求解速度慢、需要設(shè)定迭代初值等缺點(diǎn)。由于正弦函數(shù)具有無窮階的可導(dǎo)性，因而可以保證結(jié)果的收斂性。另外，若將該方法用于海底無限深的聲學(xué)半空間模型，需引入海底吸收并保證底部邊界具有足夠的深度。本文通過數(shù)值計(jì)算，驗(yàn)證了該方法在求解等聲速波導(dǎo)、聲速剖面波導(dǎo)以及分層均勻海洋波導(dǎo)中本征值與本征函數(shù)的準(zhǔn)確性，為分析均勻、非均勻海洋波導(dǎo)的聲傳播問題提供了新的思路。