劉娟 李琪

1) (哈爾濱工程大學(xué)， 水聲技術(shù)重點(diǎn)實(shí)驗(yàn)室， 哈爾濱 150001)

2) (海洋信息獲取與安全工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室(哈爾濱工程大學(xué))， 哈爾濱 150001)

3) (哈爾濱工程大學(xué)水聲工程學(xué)院， 哈爾濱 150001)

1 引 言

水平變化波導(dǎo)中的聲傳播問題是水下聲傳播研究中的熱點(diǎn)問題.受季節(jié)和日照的影響， 以及海洋中的中小尺度現(xiàn)象和內(nèi)波的作用， 海洋環(huán)境參數(shù)空間分布不均勻.作為海洋波導(dǎo)的上下界面， 海面在風(fēng)力驅(qū)動(dòng)作用下隨機(jī)起伏， 海底地形粗糙不平.聲波在水平變化波導(dǎo)中傳播會(huì)發(fā)生隨機(jī)散射， 嚴(yán)重影響能量和信息的傳輸.本文針對(duì)介質(zhì)參數(shù)和海底邊界兩種水平變化因素， 研究水平變化波導(dǎo)中的聲傳播問題.

耦合模態(tài)法是求解水平變化波導(dǎo)中聲傳播問題的常用方法之一， 由Pierce[1]和Milder[2]引入到水聲學(xué)中.Pierce-Milder 耦合模態(tài)理論將水平變化波導(dǎo)中任意位置處的聲壓表示為本地模態(tài)的疊加， 模態(tài)系數(shù)與水平位置有關(guān).通過對(duì)聲壓滿足的Helmholtz 方程在本地模態(tài)上作投影， 將Helmholtz 方程轉(zhuǎn)化為關(guān)于模態(tài)系數(shù)的二階耦合模態(tài)方程組.波導(dǎo)環(huán)境的水平變化對(duì)聲場的作用由兩個(gè)耦合矩陣描述.本地模態(tài)的計(jì)算精度直接影響耦合模態(tài)法的求解精度.對(duì)于介質(zhì)參數(shù)均勻的波導(dǎo)， 本地模態(tài)可解析計(jì)算.對(duì)于介質(zhì)參數(shù)復(fù)雜的波導(dǎo)， 需借助數(shù)值手段求解， 如譜方法[3]、多模態(tài)法[4]、有限元法[5]、有限差分法[6]及標(biāo)準(zhǔn)簡正波程序KRAKEN[7]等.Rutherford 和Hawker[8]指出， Pierce-Milder耦合模態(tài)理論在處理邊界傾斜波導(dǎo)的聲傳播問題時(shí)能量不守恒.這是由于本地模態(tài)并不符合嚴(yán)格的邊界條件.Fawcett[9]將正確的邊界條件引入到耦合模態(tài)方程組的推導(dǎo)中， 提出一種能量守恒的耦合模態(tài)理論.但該理論不僅要計(jì)算兩個(gè)耦合矩陣， 還需計(jì)算兩個(gè)界面矩陣.這些矩陣與本地模態(tài)在水平方向上的導(dǎo)數(shù)有關(guān)， 只能近似求解.另外， 直接積分求解二階耦合模態(tài)方程組可能會(huì)遇到數(shù)值發(fā)散的問題[4].本地模態(tài)和其水平導(dǎo)數(shù)的計(jì)算復(fù)雜度及二階耦合模態(tài)方程的數(shù)值求解問題使得Fawcett的方法在實(shí)際中難以廣泛應(yīng)用.Abawi[10]通過忽略高階耦合項(xiàng)和后向散射場的能量， 結(jié)合拋物方程方法推導(dǎo)得到耦合簡正波-拋物方程理論， 有效解決了耦合模態(tài)方程的求解問題.在此基礎(chǔ)上， 彭朝暉和張仁和[11]采用廣義相積分理論和波束位移射線簡正波理論計(jì)算本地模態(tài)， 實(shí)現(xiàn)了聲場的快速計(jì)算.莫亞梟等[12]通過忽略高階耦合項(xiàng)并只取前向場， 提出了一種大步長格式的二維耦合模態(tài)方法.這些單向耦合模態(tài)理論對(duì)于環(huán)境特性水平不變或緩變波導(dǎo)是正確的， 但不適用于后向場能量不可忽視的情況[13].針對(duì)求解二階耦合模態(tài)方程中可能出現(xiàn)的數(shù)值發(fā)散問題， Pagneux 等[4，14]提出一種數(shù)值收斂的多模態(tài)導(dǎo)納法， 分析了變截面聲管中的聲傳播問題.多模態(tài)導(dǎo)納法通過引入導(dǎo)納矩陣將Helmholtz 方程的邊值問題等價(jià)成兩個(gè)一階演化方程的初值問題， 采用Runge-Kutta 法[4，15]或Magnus 法[14，16，17]數(shù)值計(jì)算獲得穩(wěn)定聲場解.目前研究海底地形變化的模態(tài)耦合理論已相對(duì)成熟， 而介質(zhì)參數(shù)水平變化波導(dǎo)中的聲傳播問題仍有待研究.對(duì)于介質(zhì)參數(shù)水平緩變(相比聲波波長)的海洋環(huán)境， 絕熱近似耦合模態(tài)理論[18]通過忽略模態(tài)間的耦合作用， 得到形式簡單的耦合模態(tài)方程.然而， 對(duì)于垂直于內(nèi)波波峰方向上的聲傳播問題， 環(huán)境特性會(huì)迅速變化， 絕熱近似理論通常無效[19].因此， 有必要推導(dǎo)一種雙向耦合模態(tài)方法， 考慮介質(zhì)參數(shù)水平變化對(duì)模態(tài)耦合及聲場的影響.

前文所述均為連續(xù)耦合模態(tài)理論(continuous coupled mode method)， Evans[20]將水平變化區(qū)域分為若干段水平不變波導(dǎo)， 提出階梯耦合模態(tài)理論(stepwise coupled mode method).駱文于等[7，21]引入全局矩陣算法， 推導(dǎo)出一種數(shù)值穩(wěn)定的階梯耦合模態(tài)方法.階梯耦合模態(tài)方法數(shù)值實(shí)現(xiàn)相對(duì)簡單， 在水聲領(lǐng)域有廣泛應(yīng)用.相比階梯耦合模態(tài)理論， 連續(xù)耦合模態(tài)方法有兩個(gè)優(yōu)勢: 1)給出了耦合矩陣的具體表達(dá)式， 直觀地體現(xiàn)水平變化因素對(duì)模態(tài)間能量交換及聲場的影響; 2)方法可以延伸應(yīng)用于更實(shí)際的海洋波導(dǎo)中的聲傳播問題， 如: 三維模型[22]、含散射體[23，24]等復(fù)雜波導(dǎo)環(huán)境.

鑒于上述問題， 本文將結(jié)合Fawcett 耦合模態(tài)理論和多模態(tài)導(dǎo)納法， 提出一種能量守恒且數(shù)值收斂的耦合模態(tài)方法， 研究介質(zhì)參數(shù)和邊界水平變化波導(dǎo)中的聲傳播問題.全文結(jié)構(gòu)如下: 首先， 針對(duì)介質(zhì)參數(shù)及邊界水平變化波導(dǎo)， 推導(dǎo)基于多模態(tài)導(dǎo)納法的耦合模態(tài)方法， 給出耦合模態(tài)方程及求解耦合模態(tài)方程的Magnus 數(shù)值積分方法; 其次， 使用該方法計(jì)算水平變化波導(dǎo)中的聲場， 與COMSOL的計(jì)算結(jié)果比較， 驗(yàn)證方法的準(zhǔn)確性; 最后， 給出討論和總結(jié).

2 理論方法

針對(duì)水平變化波導(dǎo)， 本節(jié)給出基于多模態(tài)導(dǎo)納法的耦合模態(tài)方法(coupled mode method)， 下文中簡稱CMM.

2.1 雙向耦合模態(tài)理論

考慮水平方向半無窮、垂直方向受限的二維平面波導(dǎo)， 如圖1 所示.上邊界為水平不變的空氣-海水絕對(duì)軟邊界， 下邊界為水平變化的海水-巖石絕對(duì)硬邊界， 以連續(xù)函數(shù) H (x) 表 示.聲速 c (x，y) 和密度 ρ (x，y) 均 是空 間 坐標(biāo) 的連 續(xù) 函 數(shù).假 設(shè) 區(qū) 域Ω1:x ∈[0，xr]為水平變化區(qū)域， 聲波在該區(qū)域中發(fā)生散射.區(qū)域 Ω2:x ＞xr為水平不變區(qū)域， 聲波在該區(qū)域中全部透射.入射波從 x =0 處向右入射.

圖1 水平變化波導(dǎo)示意圖Fig.1.Configuration of a waveguide with range-dependent environments.

省略時(shí)間因子 e-iωt， 聲壓滿足Helmholtz 方程

和邊界條件

其中p 是聲壓， ω 是角頻率， n 表示外法線方向.

根據(jù)耦合模態(tài)理論， 波導(dǎo)中任意位置處的聲壓可用本地模態(tài)疊加求和表示， 這里本地模態(tài)是指上下邊界分別滿足絕對(duì)軟和絕對(duì)硬邊界條件， 聲速和密度分別等于本地聲速 c (y;x) 和本地密度 ρ (y;x) 的水平不變波導(dǎo)中的簡正波.而依據(jù)多模態(tài)理論， 本地模態(tài)可用一組正交完備的本地本征函數(shù)疊加求和表示， 因此， 波導(dǎo)中的聲壓可表示為

其中: Pn(x) 為 模態(tài)系數(shù); φn(y;x) 為本地本征函數(shù)，代表上下邊界分別滿足絕對(duì)軟和絕對(duì)硬邊界條件.波導(dǎo)深度為本地深度 H (x) 的均勻(聲速和密度均為常數(shù))波導(dǎo)中的簡正波.用本地本征函數(shù)φn(y;x)作基函數(shù)的優(yōu)勢在于 φn(y;x) 及其導(dǎo)數(shù)的解析表達(dá)式容易計(jì)算， 而真正的本地模態(tài)及其導(dǎo)數(shù)往往不易解析求解.φn(y;x) 的表達(dá)式為φn(y;x)=滿足正交性φn(y;x)dy =δmn.注意(3)式中的求和號(hào)表示在數(shù)值實(shí)現(xiàn)中要對(duì)級(jí)數(shù)求和作截?cái)嗵幚砥渲蠳 表示截?cái)鄶?shù).對(duì)于水平不變波導(dǎo)， 求和項(xiàng)通常只需包含全部傳播模態(tài).對(duì)于水平變化波導(dǎo)， 由于衰逝模態(tài)對(duì)水平變化區(qū)域中近場聲壓的貢獻(xiàn)不可忽略， 求和項(xiàng)中需包含全部的傳播模態(tài)和部分衰逝模態(tài)[14].

代入聲壓展開(3)式， 得到二階耦合模態(tài)方程組

寫成矩陣形式

其中， 列向量P 中元素為 Pn(x) ， 矩陣A， B， C，D， E， F， G， J 和K 中的元素分別為

(6)式中的各個(gè)矩陣描述了水平變化因素對(duì)聲場的影響.矩陣A， C 和K 是由于基函數(shù)不是真正的本地模態(tài)產(chǎn)生的系數(shù)矩陣.B， D， J 和G 是耦合矩陣， 描述模態(tài)間的耦合作用.E 和F 是邊界矩陣， 來自于計(jì)算的過程中加入的嚴(yán)格下邊界條件， 即 ?yp(H)=H′?xp(H) ， 目的是修正本地本征函數(shù)的邊界條件， 使二階耦合模態(tài)方程(6)滿足能量守恒， 詳細(xì)過程可參考文獻(xiàn)[9].

至此， Helmholtz 方程引導(dǎo)的聲傳播問題被轉(zhuǎn)化為關(guān)于模態(tài)系數(shù)P 的二階耦合模態(tài)方程的邊值問題.由于邊界矩陣E 和F 修正了邊界條件， (6)式已滿足邊界條件， 對(duì)(6)式的求解只需考慮聲源條件和輻射條件.但根據(jù)輻射條件從右向左積分求解(6)式時(shí)， 由于衰逝模態(tài)在這個(gè)方向上是指數(shù)發(fā)散的， 導(dǎo)致數(shù)值計(jì)算結(jié)果不收斂， 故而(6)式并不是理想形式.根據(jù)多模態(tài)導(dǎo)納法的思想， 將二階微分方程形式下的耦合模態(tài)方程重構(gòu)為兩個(gè)一階演化方程.定義列向量s 為

(7)式對(duì)x 作全微分， 有s′=AP′′+(A′+D)P′+D′P ， 矩陣 A′中的元素為

即A′=E+B+D+DT.矩陣 D′中的元素為

即 D′=F +G+J +L ， 其中L 中的元素為Lmn=因此， 有

根據(jù)(6)—(8)式， 有

其中用到了等式 L =DTA-1D.下面給出該等式的證明， 利用本地本征函數(shù) φn的完備性， 將函數(shù)代 入 矩 陣D 中 有寫成矩陣形式，得到 D =AWT， 其中 W 中的元素為 wnn′， 同理可得即L=W AWT， 因 此 DTA-1D =W ATA-1AWT，由于矩陣A 為對(duì)稱陣， 有DTA-1D =L.

最終， 得到兩個(gè)一階耦合模態(tài)方程

(10)式可通過引入導(dǎo)納矩陣獲得穩(wěn)定的聲場解， 并且無需求解邊界矩陣E 和F.導(dǎo)納矩陣Y 代表模態(tài)域的Dirichlet to Neumann (DtN)算子[17，25]， 定義為 s =Y P.將導(dǎo)納矩陣代入(10)式， 得到導(dǎo)納矩陣滿足的Riccati 方程

和聲壓的模態(tài)系數(shù)滿足的一階微分方程

(11)式的初值條件為區(qū)域 Ω2中的輻射條件Y(xr)=事實(shí)上， 矩陣A-1(xr)Y(xr)的本征值等于 i kxn， 其中 kxn對(duì)應(yīng)波導(dǎo)在 Ω2區(qū)域中水平方向上的波數(shù).(12)式的初值條件為 x =0 處的聲壓對(duì)應(yīng)的模態(tài)系數(shù).工程應(yīng)用及數(shù)值實(shí)現(xiàn)中，聲源通常為入射波， 而 x =0 處的聲壓包含入射波成分及反射波成分.為了計(jì)算反射波， 定義模態(tài)域的反射矩陣R 為 P-(0)=RP+(0) ， 其中 P+和P-代表前向傳播和后向傳播的模態(tài)系數(shù)分量， 有

根據(jù)(12)—(14)式， 得到

其中 P+(0) 為入射波的模態(tài)系數(shù).本文中的導(dǎo)納矩陣和聲壓模態(tài)系數(shù)均采用Magnus 數(shù)值積分算法求解， 詳細(xì)步驟見2.3 節(jié).

另外， 散射矩陣是描述散射區(qū)域性質(zhì)的主要參量， 是連接入射信息和出射信息的重要橋梁.針對(duì)圖1 所示模型， 散射矩陣包含反射矩陣R 和透射矩陣T.透射矩陣定義為 P+(xr)=T P+(0).引入傳播矩陣Q， 滿足 P (xr)=Q(x)P(x) ， 根據(jù)(10)式，有 Q′=-Q(-A-1D+A-1Y) ， 初值條件為Q(xr)等于單位陣.透射矩陣則可寫為 T =Q(0)(I +R).

總而言之， 對(duì)于一個(gè)水平變化波導(dǎo)， 雙向耦合模態(tài)方法(two-way CMM)通過輻射條件確定任意水平位置處的導(dǎo)納矩陣， 進(jìn)而計(jì)算散射區(qū)域?qū)?yīng)的反射矩陣， 再根據(jù)入射條件推出聲壓的模態(tài)展開系數(shù)， 最后代入聲壓的展開式中得到波導(dǎo)中的穩(wěn)態(tài)聲場.

2.2 單向近似與絕熱近似耦合模態(tài)理論

當(dāng)波導(dǎo)環(huán)境參數(shù)水平連續(xù)緩變(相比波長)時(shí)，單向傳播模型及絕熱耦合模態(tài)理論能夠近似求解聲場[13].本節(jié)將給出基于多模態(tài)導(dǎo)納法的單向近似(one-way CMM)和絕熱近似耦合模態(tài)(adiabatic CMM)方法.

單向近似理論假設(shè)后向散射場的能量相比前向場的能量是可忽略不計(jì)的高階小量， 等價(jià)于反射矩陣近似為零 R ≈0， 即 Y0≈Y (0) ， 表示任意位置處的導(dǎo)納矩陣無需迭代計(jì)算， 近似為Y(x)=根據(jù)(10)式， 有

初值條件為入射波對(duì)應(yīng)的模態(tài)系數(shù) P (0)=P+(0).該單向耦合模態(tài)方法避免了從無窮遠(yuǎn)輻射條件迭代計(jì)算導(dǎo)納矩陣的過程， 無需計(jì)算反射矩陣， 相比雙向耦合模態(tài)方法提高了計(jì)算速度.

絕熱近似理論假設(shè)模態(tài)間絕熱耦合， 不發(fā)生能量交換.因此， 忽略耦合矩陣項(xiàng)， 得到水平變化波導(dǎo)中的絕熱近似耦合模態(tài)方程

同樣地， 初值條件為入射波對(duì)應(yīng)的模態(tài)系數(shù)P (0)=P+(0).

2.3 Magnus 數(shù)值積分方法

本文通過Magnus 數(shù)值積分方法求解導(dǎo)納矩陣及聲壓模態(tài)系數(shù).Magnus 法是一種高效的數(shù)值積分算法， 即使是劇烈振蕩的解函數(shù)也只需稀疏的離散點(diǎn)， 這種大步長的計(jì)算特性使得Magnus 法在高頻條件下能夠?qū)崿F(xiàn)快速計(jì)算[14].Magnus 法計(jì)算導(dǎo)納矩陣與聲壓模態(tài)系數(shù)的具體步驟如下.

首先將散射區(qū)域 Ω1中的橫坐標(biāo)從右向左離散為， 其 中，=0 ， 離散間隔為＜0.將(10)式重寫為

(19)式的Magnus 數(shù)值解為

其中Magnus 矩陣Γj的形式取決于不同階的Magnus 積分方法， 二階Magnus 矩陣為

四階Magnus 矩陣為

將 eΓj寫為分塊矩陣， 即

最終得到導(dǎo)納矩陣和聲壓的模態(tài)系數(shù)

3 數(shù)值結(jié)果

本節(jié)將利用第2 節(jié)提出的耦合模態(tài)方法， 求解水平變化波導(dǎo)中分布源及點(diǎn)源聲傳播問題， 并與COMSOL 參考解比較， 驗(yàn)證耦合模態(tài)方法的正確性.針對(duì)水平緩變波導(dǎo)的聲傳播問題， 討論單向近似及絕熱近似耦合模態(tài)方法.數(shù)值驗(yàn)證雙向耦合模態(tài)方法的能量守恒特性.最后針對(duì)淺海孤立子內(nèi)波中的低頻聲傳播問題， 分析模態(tài)間的耦合作用.

3.1 節(jié)至3.3 節(jié)中考慮如圖1 所示的含聲速和密度 剖 面 的 水平 變化 波 導(dǎo)， 密度 ρ (x，y) (單 位 為kg/m3)和聲速 c (x，y) (單位為 m /s )均是空間坐標(biāo)的函數(shù)， 表達(dá)式分別為

其中， xr為水平變化距離， ρ1和 ρ2分別代表密度垂直和水平變化率， c1和 c2分別代表聲速垂直和水平變化率.上邊界水平不變， 下邊界的幾何參數(shù)為

其中 h1代表海底傾斜率.

3.1 分布源和點(diǎn)源聲傳播問題

首先計(jì)算分布源激發(fā)聲場.環(huán)境參數(shù)選取為ρ1=0.2 ， ρ2=0.2 ， c1=0.2 ， c2=0.4 ， h1=-0.1.H =200 m， xr=1600 m.聲源為從 x =0 處向右入射的分布源 pi=εφ1(y;0)=sin(1.5πy/200) ， 其中ε=10 為入射波歸一化系數(shù).頻率 f =20 Hz.在該頻率下， 波導(dǎo)在 x =0 m 處有五階可傳播模態(tài)， 在水平不變區(qū)域中有三階可傳播模態(tài).圖2(a)為利用雙向CMM((10)式)計(jì)算得到的聲壓幅值(Pa)分布， 截?cái)鄶?shù)選取為 N =10 ， 圖2(b)為使用COMSOL 計(jì)算得到的聲壓幅值(Pa)分布， 圖2(c)為深度 y =20 m 處聲壓幅值沿x 軸的分布， 圖2(d)為深度 y =50 m 處聲壓幅值沿x 軸的分布.可以看出水平變化區(qū)域中存在明顯的后向散射效應(yīng)， 模態(tài)間耦合作用劇烈， 兩種方法的計(jì)算結(jié)果一致.

圖2 水平變化波導(dǎo)中的聲場(聲源為分布源， 頻率為20 Hz) (a) 利用CMM 計(jì)算得到的聲壓幅值分布， 截?cái)鄶?shù) N =10 ; (b) 利用COMSOL 計(jì)算得到的聲壓幅值分布; (c) 深度為20 m 處， 聲壓幅值的水平分布; (d) 深度為50 m 處， 聲壓幅值的水平分布Fig.2.Sound fields in a range-dependent waveguide (the source is a distributed source at 20 Hz): (a) Sound field computed by CMM where the truncation number is N =10 ; (b) sound field computed with COMSOL; (c) sound field distribution along x direction at depth 20 m; (d) sound field distribution along x direction at depth 50 m.

CMM 方法處理點(diǎn)源(相當(dāng)于對(duì)稱三維問題中的線源)聲傳播問題的關(guān)鍵在于用本地本征函數(shù)展開表示點(diǎn)源傳播函數(shù)—Green 函數(shù)， 即

對(duì)于水平變化波導(dǎo)， 相應(yīng)的模態(tài)系數(shù)g 可寫為[24，26]

圖3 所示為點(diǎn)源激發(fā)聲場.水平變化波導(dǎo)的環(huán)境參數(shù)為 ρ1=0.2 ， ρ2=0.2 ， c1=0.2 ， c2=0.4 ，h1=-0.05.H =200 m， xr=1600 m.聲源位于(0，10) m 處.頻率為 f =20 H z.圖3(a)為利用雙向CMM((10)式)計(jì)算得到的聲壓幅值(Pa)， 截?cái)鄶?shù)選取為 N =50 ， 圖3(b)為使用COMSOL 計(jì)算得到的聲壓幅值(Pa)， 圖3(c)為深度 y =71 m 處聲壓幅值沿x 軸的分布， 圖3(d)為深度 y =101 m處聲壓幅值沿x 軸的分布.可以看出明顯的后向散射作用.兩種方法的計(jì)算結(jié)果一致， 說明雙向CMM能夠準(zhǔn)確計(jì)算水平變化波導(dǎo)中的聲場.兩種方法的計(jì)算偏差來源于兩方面: 1) CMM 在數(shù)值實(shí)現(xiàn)中對(duì)級(jí)數(shù)求和(3)式作截?cái)嗵幚恚?導(dǎo)致的部分和與真值之間的誤差; 2)本文中采用Clenshaw-Curtis 數(shù)值積分方法[26]計(jì)算系數(shù)矩陣及耦合矩陣， 離散點(diǎn)采用Chebyshev 插值點(diǎn)， 而COMSOL 計(jì)算采用三角形網(wǎng)格， CMM 與COMSOL 離散方式不同導(dǎo)致兩種結(jié)果出現(xiàn)偏差.

圖3 水平變化波導(dǎo)中的聲場(聲源為位于 ( 0，10) m 處的點(diǎn)源， 頻率為20 Hz) (a) 利用雙向CMM 計(jì)算得到的聲壓幅值分布，截?cái)鄶?shù) N =50 ; (b) 利用COMSOL 計(jì)算得到的聲壓幅值分布; (c) 深度為71 m 處， 聲壓幅值的水平分布; (d) 深度為101 m 處， 聲壓幅值的水平分布Fig.3.Sound fields in a range-dependent waveguide generated by a point source at ( 0，10) m (the frequency is 20 Hz): (a) Sound field computed by CMM where the truncation number is N =50 ; (b) sound field computed with COMSOL; (c) sound field distribution along x direction at depth 71 m; (d) sound field distribution along x direction at depth 101 m.

3.2 水平緩變波導(dǎo): 單向近似與絕熱近似

水平緩變是指環(huán)境參數(shù)在一個(gè)波長內(nèi)的水平變化量遠(yuǎn)小于波長， 對(duì)于本節(jié)的計(jì)算模型， 水平緩變代表 ρ1?1 ， c1?1 及 h1?1.考慮水平緩變波導(dǎo)算例， 假設(shè) ρ1=0.01 ， ρ2=0.2 ， c1=0.01 ， c2=0.4 ， h1=-0.001 ， H =200 m， xr=1600 m.聲源為從 x =0 處向右入射的分布源pi=εφ1(y;0)=sin(1.5πy/200).頻率 f =20 Hz.圖4(a)—(c)分別為利用雙向CMM 理論((10)式)、單向近似CMM理論((17)式)和絕熱近似CMM 理論((18)式)計(jì)算得到的聲壓幅值(Pa)分布， 圖4(d)和圖4(e)分別為接收點(diǎn)在60 m 和120 m 深度處， 聲壓幅值的水平分布.三種方法的數(shù)值實(shí)現(xiàn)采用相同的離散點(diǎn)和截?cái)鄶?shù)N = 10.可以看出， 單向近似和絕熱近似耦合模態(tài)理論均可較為準(zhǔn)確地計(jì)算水平緩變波導(dǎo)中的近距離聲場， 但計(jì)算誤差隨水平距離的遞推逐漸累積.

3.3 能量守恒驗(yàn)證

本節(jié)討論數(shù)值聲場解的能流守恒， 理論驗(yàn)證將在第4 節(jié)討論部分中給出.根據(jù)能流定義

代入聲壓展開式(3)， 寫成矩陣形式， 得到模態(tài)(本地本征函數(shù))域的能流表達(dá)式

其中上標(biāo)H 代表共軛轉(zhuǎn)置.

圖4 水平緩變波導(dǎo)中的聲場(聲源為分布源， 頻率為20 Hz) (a) 利用雙CMM 計(jì)算得到的聲壓幅值分布; (b) 利用單向近似CMM 計(jì)算得到的聲壓幅值分布; (c) 利用絕熱近似CMM 計(jì)算得到的聲壓幅值分布; (d) 深度為60 m 處， 聲壓幅值的水平分布;(e) 深度為120 m 處， 聲壓幅值的水平分布Fig.4.Sound fields in a waveguide with weak range dependence generated by a distributed source (the frequency is 20 Hz): (a) Sound field computed by two-way CMM; (b) sound field computed with one-way CMM; (c) sound field computed with adiabatic CMM;(d) sound field distribution along x direction at depth 60 m; (e) sound field distribution along x direction at depth 120 m.

為直觀起見， 考慮邊界水平不變， 密度均勻，聲速水平變化波導(dǎo)中的聲傳播.假設(shè) ρ1=0 ，ρ2=0 ， c1=0.1 ， c2=0 ， h1=0 ， H =200 m.頻 率f =10 Hz.在 x =0 m 處波導(dǎo)中有三階可傳播模態(tài).隨著聲速的遞增， 在水平不變區(qū)域中僅有兩階可傳播模態(tài).由于介質(zhì)參數(shù)與y 無關(guān)， 本地本征函數(shù)即為實(shí)際本地簡正波.聲源為從 x =0 m 處向右入射的分布源 pi=εφ1(y;0)=sin(1.5πy/200) ， 即第二階本地簡正波.圖5(a)為雙向CMM ((10)式)計(jì)算得到的聲壓幅值(Pa)分布， 圖5(b)為聲壓的模態(tài)系數(shù) | Pn(x)| 的分布， 圖5(c)為歸一化能流E(x)/E(0).從圖5(c)可以看出， 能流為常數(shù)， 代表數(shù)值解符合能量守恒.由于衰逝模態(tài)不傳播能量，圖5(b)中僅畫出可傳播模態(tài)對(duì)應(yīng)的模態(tài)系數(shù).入射波為第二階模態(tài)， 即 P1，+(0)=1 ， 圖中顯示僅有第二階模態(tài)在波導(dǎo)中傳播， 表明聲速的水平變化沒有導(dǎo)致模態(tài)間的耦合.| P1(0)|=1.0005 ， 說明后向散射場的能量很小.從反射矩陣也可知，||RP||2～O(10-4)?1.另外， 密度均勻時(shí)， 能流E(x)∝Re聲速隨傳播距離增大， | P1(x)| 也隨之增大， 符合能量守恒定律，與數(shù)值計(jì)算結(jié)果相符.

接著在圖5 密度均勻、聲速水平變化波導(dǎo)的基礎(chǔ)上加入邊界變化.假設(shè) ρ1=0 ， ρ2=0 ， c1=0.1 ，c2=0 ， h1=-0.03.H =200 m， xr=1600 m.頻率為 10 H z.在 x =0 m 處波導(dǎo)中有三階可傳播模態(tài)， 在水平不變區(qū)域中有兩階可傳播模態(tài).聲源依舊為從 x =0 處向右入射的分布源pi=εφ1(y;0)=sin(1.5πy/200)， 即第二階本地簡正波.圖6(c)顯示數(shù)值解滿足能量守恒.圖6(b)說明雖然入射波為第二階簡正波， 但環(huán)境水平變化導(dǎo)致相鄰位置的本地簡正波之間發(fā)生耦合， 在 x =0 處前三階模態(tài)在波導(dǎo)中同時(shí)傳播， 隨著聲波的傳播， 第三階本地模態(tài)被截止， 最終在水平不變區(qū)域中僅有前兩階模態(tài)傳播.此外， | |RP||2=0.0277 ， 相比邊界不變波導(dǎo)(圖5)， 邊界水平變化較易產(chǎn)生不可忽視的后向散射作用.

3.4 淺海孤立子內(nèi)波

孤立子內(nèi)波是淺海中常見的聲學(xué)現(xiàn)象.內(nèi)波可能帶來大尺度的水體垂直位移， 導(dǎo)致低頻聲信號(hào)在傳播中產(chǎn)生波動(dòng)[27].由于每個(gè)孤立子內(nèi)波以各自的相速度傳播， 多個(gè)內(nèi)波間的相互作用復(fù)雜且難以描述內(nèi)波對(duì)模態(tài)耦合的影響.因此， 以單個(gè)孤立子內(nèi)波環(huán)境為例， 分析模態(tài)間的耦合作用.

考慮圖7 所示的包含單個(gè)孤立子內(nèi)波的淺海波導(dǎo)， 海面為絕對(duì)軟界面， 海底為絕對(duì)硬界面.水體深度為60 m， 依據(jù)聲速分布被分為三層， 上層為等聲速層， 聲速為 c1=1530 m/s， 中間為內(nèi)波層，聲速為

下層為等聲速層， 聲速為 c1=1490 m/s.各層中密度均勻 ρ =1000 kg/m3.內(nèi)波的幾何形狀為

圖5 密度均勻、聲速水平變化波導(dǎo)中的聲傳播(聲源為第二階模態(tài)， 頻率為10 Hz) (a)聲壓幅值分布; (b)聲壓的模態(tài)系數(shù)的水平分布; (c)歸一化能流Fig.5.Sound propagation in a waveguide with constant mass density and range-dependent sound speed (the source is the second local mode， and its frequency is 10 Hz): (a) Sound field; (b) modal amplitudes distribution; (c) normalized energy flux distribution.

圖6 密度均勻、聲速及邊界水平變化波導(dǎo)中的聲傳播(聲源為第二階簡正波， 頻率為10 Hz) (a)聲壓幅值分布; (b)聲壓的模態(tài)系數(shù)的水平分布; (c)歸一化能流Fig.6.Sound propagation in a waveguide with constant mass density and range-dependent sound speed and boundary (the source is the second local mode and its frequency is 10 Hz): (a) Sound field; (b) modal amplitudes distribution; (c) normalized energy flux distribution.

圖7 孤立子內(nèi)波波導(dǎo)模型Fig.7.Configuration of a waveguide with an internal solitary wave.

這種類駝峰的內(nèi)波形狀是Korteweg-de Vries 方程的形式解， 描述弱非線性內(nèi)波的傳播[19].

研究波導(dǎo)中的模態(tài)耦合作用時(shí)， 常見且實(shí)用的手段是考慮單階簡正波激發(fā)的聲場中每個(gè)水平位置處各階簡正波的分量.雖然本文提出的耦合模態(tài)法采用的展開基函數(shù)((3)式)不是真正的本地簡正波， 展開系數(shù)也不是真正的各階簡正波的分量， 然而， 由(12)式的輻射條件可知， 矩陣的特征向量對(duì)應(yīng)各階本地簡正波在本地本征函數(shù) φn上的模態(tài)系數(shù)， 因此， 本地簡正波及聲場在本地簡正波上的分量可通過正交基變換獲得.據(jù)此， 我們數(shù)值計(jì)算得到圖7 模型中 x =0 處的真正本地簡正波， 示于圖8.

圖9(a)為入射第一階本地簡正波(即圖8(a))，頻率為100 Hz 時(shí)， 采用雙向CMM((10)式)計(jì)算得到的聲壓幅值分布(Pa).圖9(b)為聲壓在各個(gè)水平位置處各階本地簡正波的分量.顯然地， 在聲速水平不變區(qū)域中， 簡正波之間沒有耦合.在聲速水平變化區(qū)域， 簡正波之間產(chǎn)生耦合， 能量從第一階簡正波轉(zhuǎn)移到高階簡正波中.另外， 連續(xù)水平變化的聲速?zèng)]有導(dǎo)致強(qiáng)烈的反射， 故而單向傳播模型同樣適用于求解該環(huán)境中聲場.

4 討論部分

首先理論驗(yàn)證雙向CMM 滿足能量守恒.假設(shè) p (x，y) 和 p*(x，y) 為Helmholtz 方 程 的 聲 場 解，也就是

圖8 圖7 所示模型中 x =0 m 處的簡正波Fig.8.Local mode shape of the waveguide at x =0 m shown in Fig.7.

圖9 淺海孤立子內(nèi)波波導(dǎo)中的聲傳播(聲源為第一階簡正波， 頻率為100 Hz) (a)聲壓幅值分布; (b)前五階簡正波分量Fig.9.Sound propagation in a waveguide with an internal solitary wave (the incident wave is the first mode at 100 Hz): (a) Sound field; (b) modal amplitudes distribution.

(35)式 × p*- (36)式 × p ， 得到

(37)式對(duì)y 積分， 得到

代入聲壓展開式(3)， 根據(jù)(7)式， 有

因此，

(39)式對(duì)應(yīng)聲場互易性公式， (40)式對(duì)應(yīng)能流守恒公式.可見， 雙向耦合模態(tài)方法滿足能量守恒定律.

本文提出的雙向耦合模態(tài)方法為半解析半數(shù)值解法， 在推導(dǎo)過程中沒有引入任何近似假設(shè).然而， 水平變化因素—聲速、密度及邊界幾何函數(shù)均默認(rèn)為連續(xù)且可導(dǎo).對(duì)于水平變化區(qū)域中邊界導(dǎo)數(shù)不存在的情況， 可以通過保角變換[28]將復(fù)雜形狀邊界下的邊界條件轉(zhuǎn)換為水平不變邊界下的邊界條件， 將均勻Helmholtz 方程變?yōu)楹凶兓凵渎实腍elmholtz 方程， 再采用耦合模態(tài)法處理.該雙向耦合模態(tài)方法可以用于研究水平分層波導(dǎo)中的聲傳播問題， 即介質(zhì)參數(shù)不連續(xù)變化的波導(dǎo)環(huán)境， 只需恰當(dāng)?shù)丶尤虢唤缑嫣幍倪B續(xù)性條件.對(duì)于下層半無窮波導(dǎo)， 如Pekeris 波導(dǎo)， 第3 節(jié)中指出聲速水平連續(xù)緩慢變化時(shí)， 后向散射場幾乎可忽略， 因此可以通過在下層空間的聲速上加入隨深度逐漸增大的吸收項(xiàng)[29]， 近似分析深度半無窮波導(dǎo)中的聲傳播問題.

5 結(jié) 論

本文提出一種基于多模態(tài)導(dǎo)納法的耦合模態(tài)方法以研究水平變化波導(dǎo)中的聲傳播問題.根據(jù)傳統(tǒng)耦合模態(tài)理論和多模態(tài)法， 將每個(gè)位置處的聲壓用一組正交本地本征函數(shù)展開， 對(duì)Helmholtz 方程在本地本征函數(shù)上作投影， 推導(dǎo)出模態(tài)系數(shù)滿足的二階耦合模態(tài)方程組.為解決二階耦合方程組數(shù)值計(jì)算發(fā)散的問題， 引入導(dǎo)納矩陣， 將二階耦合模態(tài)方程組簡化為兩個(gè)一階演化方程， 并利用Magnus數(shù)值積分方法數(shù)值求解.利用該耦合模態(tài)方法數(shù)值求解水平變化波導(dǎo)中的聲場， 與COMSOL 參考解比較， 結(jié)果吻合， 表明該耦合模態(tài)方法能夠精確求解水平變化波導(dǎo)中的分布源及點(diǎn)源聲傳播問題.單向近似和絕熱近似耦合模態(tài)方法可以近似求解水平緩變波導(dǎo)的聲場.雙向耦合模態(tài)方法在推導(dǎo)過程中沒有任何近似， 計(jì)算誤差來源于數(shù)值實(shí)現(xiàn)， 并且Magnus 數(shù)值積分方法具有大步長的計(jì)算特性， 滿足能量守恒且數(shù)值解無條件穩(wěn)定.總之， 雙向耦合模態(tài)法方法充分考慮了波導(dǎo)環(huán)境對(duì)模態(tài)耦合的作用， 能夠精確有效求解水平變化波導(dǎo)中的聲傳播問題， 可以用于求解實(shí)際聲傳播問題.