一種水平變化波導中聲傳播問題的耦合模態(tài)法

劉娟 李琪

1) (哈爾濱工程大學， 水聲技術(shù)重點實驗室， 哈爾濱 150001)

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

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

1 引 言

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

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

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

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

2 理論方法

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

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

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

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

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

和邊界條件

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

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

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

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

寫成矩陣形式

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

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

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

(7)式對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 為對稱陣， 有DTA-1D =L.

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

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

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

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

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

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

另外， 散射矩陣是描述散射區(qū)域性質(zhì)的主要參量， 是連接入射信息和出射信息的重要橋梁.針對圖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).

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

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

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

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

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

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

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

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

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

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

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

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

四階Magnus 矩陣為

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

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

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

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

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

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

其中 h1代表海底傾斜率.

3.1 分布源和點源聲傳播問題

首先計算分布源激發(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.在該頻率下， 波導在 x =0 m 處有五階可傳播模態(tài)， 在水平不變區(qū)域中有三階可傳播模態(tài).圖2(a)為利用雙向CMM((10)式)計算得到的聲壓幅值(Pa)分布， 截斷數(shù)選取為 N =10 ， 圖2(b)為使用COMSOL 計算得到的聲壓幅值(Pa)分布， 圖2(c)為深度 y =20 m 處聲壓幅值沿x 軸的分布， 圖2(d)為深度 y =50 m 處聲壓幅值沿x 軸的分布.可以看出水平變化區(qū)域中存在明顯的后向散射效應， 模態(tài)間耦合作用劇烈， 兩種方法的計算結(jié)果一致.

圖2 水平變化波導中的聲場(聲源為分布源， 頻率為20 Hz) (a) 利用CMM 計算得到的聲壓幅值分布， 截斷數(shù) N =10 ; (b) 利用COMSOL 計算得到的聲壓幅值分布; (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 方法處理點源(相當于對稱三維問題中的線源)聲傳播問題的關(guān)鍵在于用本地本征函數(shù)展開表示點源傳播函數(shù)—Green 函數(shù)， 即

對于水平變化波導， 相應的模態(tài)系數(shù)g 可寫為[24，26]

圖3 所示為點源激發(fā)聲場.水平變化波導的環(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)式)計算得到的聲壓幅值(Pa)， 截斷數(shù)選取為 N =50 ， 圖3(b)為使用COMSOL 計算得到的聲壓幅值(Pa)， 圖3(c)為深度 y =71 m 處聲壓幅值沿x 軸的分布， 圖3(d)為深度 y =101 m處聲壓幅值沿x 軸的分布.可以看出明顯的后向散射作用.兩種方法的計算結(jié)果一致， 說明雙向CMM能夠準確計算水平變化波導中的聲場.兩種方法的計算偏差來源于兩方面: 1) CMM 在數(shù)值實現(xiàn)中對級數(shù)求和(3)式作截斷處理， 導致的部分和與真值之間的誤差; 2)本文中采用Clenshaw-Curtis 數(shù)值積分方法[26]計算系數(shù)矩陣及耦合矩陣， 離散點采用Chebyshev 插值點， 而COMSOL 計算采用三角形網(wǎng)格， CMM 與COMSOL 離散方式不同導致兩種結(jié)果出現(xiàn)偏差.

圖3 水平變化波導中的聲場(聲源為位于 ( 0，10) m 處的點源， 頻率為20 Hz) (a) 利用雙向CMM 計算得到的聲壓幅值分布，截斷數(shù) N =50 ; (b) 利用COMSOL 計算得到的聲壓幅值分布; (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 水平緩變波導: 單向近似與絕熱近似

水平緩變是指環(huán)境參數(shù)在一個波長內(nèi)的水平變化量遠小于波長， 對于本節(jié)的計算模型， 水平緩變代表 ρ1?1 ， c1?1 及 h1?1.考慮水平緩變波導算例， 假設(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)式)計算得到的聲壓幅值(Pa)分布， 圖4(d)和圖4(e)分別為接收點在60 m 和120 m 深度處， 聲壓幅值的水平分布.三種方法的數(shù)值實現(xiàn)采用相同的離散點和截斷數(shù)N = 10.可以看出， 單向近似和絕熱近似耦合模態(tài)理論均可較為準確地計算水平緩變波導中的近距離聲場， 但計算誤差隨水平距離的遞推逐漸累積.

3.3 能量守恒驗證

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

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

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

圖4 水平緩變波導中的聲場(聲源為分布源， 頻率為20 Hz) (a) 利用雙CMM 計算得到的聲壓幅值分布; (b) 利用單向近似CMM 計算得到的聲壓幅值分布; (c) 利用絕熱近似CMM 計算得到的聲壓幅值分布; (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.

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

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

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

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

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

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

圖5 密度均勻、聲速水平變化波導中的聲傳播(聲源為第二階模態(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 密度均勻、聲速及邊界水平變化波導中的聲傳播(聲源為第二階簡正波， 頻率為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)波波導模型Fig.7.Configuration of a waveguide with an internal solitary wave.

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

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

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

4 討論部分

首先理論驗證雙向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)波波導中的聲傳播(聲源為第一階簡正波， 頻率為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)式對y 積分， 得到

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

因此，

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

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

5 結(jié) 論

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