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

劉娟 李琪

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

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

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

1 引 言

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

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

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

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

2 理論方法

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

2.1 雙向耦合模態理論

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

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

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

和邊界條件

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

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

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

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

寫成矩陣形式

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

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

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

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

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

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

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

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

最終， 得到兩個一階耦合模態方程

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

和聲壓的模態系數滿足的一階微分方程

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

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

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

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

總而言之， 對于一個水平變化波導， 雙向耦合模態方法(two-way CMM)通過輻射條件確定任意水平位置處的導納矩陣， 進而計算散射區域對應的反射矩陣， 再根據入射條件推出聲壓的模態展開系數， 最后代入聲壓的展開式中得到波導中的穩態聲場.

2.2 單向近似與絕熱近似耦合模態理論

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

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

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

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

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

2.3 Magnus 數值積分方法

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

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

(19)式的Magnus 數值解為

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

四階Magnus 矩陣為

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

最終得到導納矩陣和聲壓的模態系數

3 數值結果

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

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

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

其中 h1代表海底傾斜率.

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

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

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

對于水平變化波導， 相應的模態系數g 可寫為[24，26]

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

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

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

3.3 能量守恒驗證

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

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

其中上標H 代表共軛轉置.

圖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.

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

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

3.4 淺海孤立子內波

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

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

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

圖5 密度均勻、聲速水平變化波導中的聲傳播(聲源為第二階模態， 頻率為10 Hz) (a)聲壓幅值分布; (b)聲壓的模態系數的水平分布; (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)聲壓的模態系數的水平分布; (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 孤立子內波波導模型Fig.7.Configuration of a waveguide with an internal solitary wave.

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

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

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

4 討論部分

首先理論驗證雙向CMM 滿足能量守恒.假設 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 淺海孤立子內波波導中的聲傳播(聲源為第一階簡正波， 頻率為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)， 根據(7)式， 有

因此，

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

本文提出的雙向耦合模態方法為半解析半數值解法， 在推導過程中沒有引入任何近似假設.然而， 水平變化因素—聲速、密度及邊界幾何函數均默認為連續且可導.對于水平變化區域中邊界導數不存在的情況， 可以通過保角變換[28]將復雜形狀邊界下的邊界條件轉換為水平不變邊界下的邊界條件， 將均勻Helmholtz 方程變為含有變化折射率的Helmholtz 方程， 再采用耦合模態法處理.該雙向耦合模態方法可以用于研究水平分層波導中的聲傳播問題， 即介質參數不連續變化的波導環境， 只需恰當地加入交界面處的連續性條件.對于下層半無窮波導， 如Pekeris 波導， 第3 節中指出聲速水平連續緩慢變化時， 后向散射場幾乎可忽略， 因此可以通過在下層空間的聲速上加入隨深度逐漸增大的吸收項[29]， 近似分析深度半無窮波導中的聲傳播問題.

5 結 論

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