潘俊杰，劉煥文，李長江

（1.浙江海洋大學(xué) 船舶與海運學(xué)院，浙江 舟山 316022；2.浙江微昊工程技術(shù)有限公司，浙江 舟山 316000）

引 言

Bragg 共振(Bragg resonance)是一類特殊的共振現(xiàn)象，最早由Bragg 父子于1913 年在對X 射線通過晶體引起的反射進行實驗時所發(fā)現(xiàn)[1].通過理論分析，他們建立了著名的Bragg 共振原理：當(dāng)晶體間距為X 射線波長一半的整數(shù)倍時，X 射線的反射達到最大.這一原理在固體力學(xué)甚至整個物理學(xué)中得到廣泛應(yīng)用，如光纖光柵傳感器和光纖光柵濾波器等，而Bragg 父子也因為此重要發(fā)現(xiàn)于1915 年雙雙獲得諾貝爾物理學(xué)獎.

Bragg 共振現(xiàn)象后來被發(fā)現(xiàn)在其他很多領(lǐng)域同樣存在.至于在水波領(lǐng)域的拓廣，主要始于Davies[2]的工作.他利用常規(guī)攝動分析，導(dǎo)出了水波越過正弦沙紋海底時Bragg 共振的激發(fā)條件：正弦沙紋的波長為表面入射波波長的一半.該結(jié)果同時得到了實驗數(shù)據(jù)[3-4]的印證.這些早期的開創(chuàng)性工作引發(fā)了大量有關(guān)正弦沙紋地形激發(fā)Bragg 共振反射的后繼性研究[5-19].

受正弦沙紋海底激發(fā)Bragg 共振形成駐波護岸甚至駐波造陸的啟示，Mei 等[20]于1988 年首次提出人工沙壩的概念，建議在近岸波浪破碎帶之外修建與岸平行且等間距周期排列的人工沙壩.相比于傳統(tǒng)防波堤，系列人工沙壩具有諸多優(yōu)點：第一，因沒于水下，沙壩兩側(cè)的水體和魚類等生物均可以不受阻礙地自由交換，而且堤上容易長出苔蘚，可吸引浮游生物的聚集，有利于改善海域環(huán)境和保持海洋生態(tài)系統(tǒng).第二，沙壩浸沒水下，不會產(chǎn)生視覺沖擊，破壞近海景觀環(huán)境.第三，壩頂可設(shè)計成低于平均低潮位，小型船舶航行不受影響.第四，多個小沙壩散布在較寬海床，不會給不夠堅硬的海底增加大的承重.最后，科學(xué)合理地布置沙壩，使之滿足Bragg 共振條件，即沙壩間距約為當(dāng)?shù)厝肷洳úㄩL一半，則能對入射波產(chǎn)生強反射.

人工沙壩的概念自提出以來，受到了學(xué)者們的廣泛關(guān)注[21-26]，2015 年前的研究進展可參見文獻[27].2016 年，Liu 等[28]除了采用物理實驗研究半圓形人工沙壩的Bragg 共振，同時針對半圓形的特殊形狀運用多極展開和極坐標(biāo)變換，使得對半圓形沙壩上方的變水深地形不需采用分段常水深的臺階系列去近似逼近，也可以直接進行特征函數(shù)展開.2017 年，Zeng 等[29]建立了線性長波越過斜坡上系列矩形沙壩時反射系數(shù)的解析解.2020 年，Liu 等[30]構(gòu)造了線性水波越過梯形人工沙壩時修正緩坡方程的Taylor 級數(shù)形式的解析解.他們第一次揭示了當(dāng)梯形沙壩相對于入射波長的無維化上底寬、下底寬以及無維化水深和浸沒深固定時，反射系數(shù)是關(guān)于沙壩間距相對于表面入射波波長的無量綱變量兩倍的周期函數(shù)，且周期為1.這說明海底地形波動和表面水波波動的周期性能夠通過反射系數(shù)的周期性表現(xiàn)出來.他們還發(fā)現(xiàn)Bragg 主振的振幅并非總是大于其他所有高階Bragg 共振的振幅.另外，隨著梯形沙壩寬度或高度的增加，Bragg 主振峰和高階共振峰的相位向低頻下移的現(xiàn)象變得越來越明顯.

最近, 受拋物形鏡面對太陽光線聚焦作用的啟發(fā), Zhang 等[31]首次將直線平行的周期波動地形改進設(shè)計成V-形彎曲后的周期波動地形, 使之可同時形成水波的Bragg 共振和聚焦雙重效應(yīng), 為增大波能的提取效率開啟了新的思路.Hao 等[32]在此基礎(chǔ)上將其拓廣設(shè)計成拋物形彎曲后的周期波動地形.

海底地形中除了周期排列的正弦沙紋和周期排列的系列人工沙壩可引起水波的Bragg 共振反射以外，互相近似平行且周期排列的溝槽系列同樣會引起水波的Bragg 共振反射.最近，有關(guān)水波越過周期排列的系列溝槽所導(dǎo)致的Bragg 共振反射已有學(xué)者開始進行研究.2018 年，Kar 等[33]采用邊界元法數(shù)值研究了水波越過海底溝槽的情形.2020 年，Kar 等[34]研究了線性淺水波越過海底斜坡上系列矩形溝槽時的Bragg 共振，結(jié)果顯示，帶有海底斜坡的系列矩形溝槽比不帶海底斜坡的系列矩形溝槽導(dǎo)致的Bragg 反射要強.2021 年，Guo等[35]基于Liu 等[36]建立的波浪被單個旋輪線形潛堤和溝槽反射的修正緩坡方程解析解，進一步解析研究了旋輪線形Bragg 潛堤和溝槽對線性波的Bragg 共振反射，第一次發(fā)現(xiàn)了當(dāng)水波越過周期排列的一列溝槽時，Bragg 共振反射的主振相位及高階共振相位均出現(xiàn)頻率上移的現(xiàn)象.這與正弦沙紋地形或周期人工沙壩激發(fā)的Bragg 共振所呈現(xiàn)的特點截然不同.在此之前我們熟知的是，當(dāng)水波越過海底正弦沙紋地形或周期排列的系列沙壩時，Bragg 共振反射的主振及高階共振的激發(fā)頻率均小于傳統(tǒng)Bragg 原理所預(yù)測的頻率，稱之為頻率下移(frequency downshift)或相位下移(phase downshift).頻率下移現(xiàn)象引起了眾多學(xué)者的興趣和關(guān)注，對正弦沙紋地形激發(fā)的頻率下移現(xiàn)象的研究參見文獻[8-9, 11]，對人工沙壩激發(fā)的頻率下移現(xiàn)象的研究可參見文獻[24-25, 28, 30, 37].最近，Liu 等[38]基于Bloch 波的能帶理論研究了線性淺水波越過拋物形和半余弦形兩類人工沙壩的頻率下移，而Liang 等[18]則基于Mathieu 不穩(wěn)定定理給出了預(yù)測主振偏移性質(zhì)的定量公式.Guo等[35]所發(fā)現(xiàn)的與以往完全不同的頻率上移現(xiàn)象表明，周期排列的系列溝槽和系列沙壩分別激發(fā)表面波產(chǎn)生Bragg 共振反射的激發(fā)頻率有著本質(zhì)的不同，值得進一步研究.

由于Guo 等[35]僅僅研究了水波被周期排列的有限個旋輪線形溝槽反射的情形，那么一個自然而然的問題就是，Bragg 共振反射峰值的相位上移現(xiàn)象是否對有限周期排列的任意形狀的溝槽而言都是普遍現(xiàn)象呢？換言之，到底是有限周期排列的任意形狀溝槽都將引起相位上移現(xiàn)象，還是僅僅只有具旋輪線形狀的有限周期系列溝槽才會引起相位上移現(xiàn)象？為了回答此問題，本文中我們基于修正緩坡方程[10]，擬進一步解析研究線性波越過海底周期排列的有限個拋物形溝槽時的Bragg 共振反射.首先，我們將每個拋物形溝槽所在的區(qū)間分為兩個單調(diào)子區(qū)間.然后在每個小區(qū)間上引入變量替換，使得帶隱函數(shù)系數(shù)的修正緩坡方程轉(zhuǎn)變?yōu)橄禂?shù)全部為顯函數(shù)的等價微分方程，據(jù)此構(gòu)造了等價微分方程的Frobenius 級數(shù)解，并確定了級數(shù)解的收斂區(qū)域.最后，利用質(zhì)量守恒耦合條件和兩個相鄰子區(qū)間的公共邊界的耦合，建立了N個周期排列的拋物形溝槽引起表面波反射的反射系數(shù)的解析公式.基于解析公式，討論了拋物形溝槽的個數(shù)、深度和寬度對Bragg 共振反射的影響.

1 拋物形系列溝槽引起線性水波反射的解析解

考慮海底地形為如圖1 所示的N個周期排列的全等拋物形溝槽，其水深函數(shù)如下：

圖1 有限個周期排列的全等拋物形溝槽示意圖Fig.1 A schematic diagram for a finite periodic array of identical parabolic trenches

我們研究一列線性水波在N個全等溝槽上的傳播.假定拋物形溝槽較淺，底部變化較緩，比如最大斜率不超過1∶1，則根據(jù)文獻[10]，自由水面位移函數(shù)η(x)滿足如下修正緩坡方程：

其中

k(x)為波數(shù)，由如下隱式的線性波色散關(guān)系決定：

其中 ω為角頻率，g為重力加速度.

不妨設(shè)具有單位振幅的波浪自左向右傳播.易知，在所有常水深區(qū)域，修正緩坡方程(2)退化為Helmholtz方程，因此相應(yīng)的自由液面位移函數(shù)容易解出為

其中k0為水深h0對應(yīng)的波數(shù)，AR和AT分別為反射和透射系數(shù)，所有復(fù)常數(shù)AR，AT，A(5j)，A(6j)均待定.

接下來致力于求解變水深溝槽區(qū)域上的自由液面位移函數(shù).在第j個溝槽的前壁區(qū)間Ij,1=[xj,1,xj,2]上，注意水深h(x)為單調(diào)函數(shù)，所以存在反函數(shù)：

運用Liu 和Zhou[39]與Liu 和Xie[40]建立的求解修正緩坡方程的變量替換技術(shù)，在修正緩坡方程(2)中引入如下變量替換：

則利用線性波色散關(guān)系(6)的等價形式

可將式(8)改寫為

注意到水深函數(shù)(1)在區(qū)間Ij,1上嚴(yán)格單調(diào)遞減，易證，由方程(9)和(11)定義的變量替換x→K為1-1 映射.顯然，該映射將Ij,1映射到區(qū)間[K0,K1]=[k0h0,k1h1]，其中k1為水深h1對應(yīng)的波數(shù).原來的位移函數(shù)η(x)在變量替換x→K后變?yōu)樾伦兞縆的函數(shù)，記為ξ(K)，即

將式(8)代入文獻[39]中的顯式方程(19)，則本文中關(guān)于舊變量x的隱式修正緩坡方程(2)被成功轉(zhuǎn)化為如下關(guān)于新變量K的顯式常微分方程：

其中

式中

方程(13)盡管已經(jīng)由隱式常微分方程轉(zhuǎn)變?yōu)轱@式常微分方程，但其系數(shù)均為變系數(shù)，非常冗長復(fù)雜，幾乎不可能構(gòu)造其封閉解.我們退而求其次，接下來致力于構(gòu)造它的級數(shù)解.為此需要先尋找級數(shù)展開式的展開點，它必須取在方程(13)的正則奇點，而不能取在方程(13)的本性奇點.容易看出K=K1導(dǎo)致S(K)和T(K)表達式中的分母為零，顯然為方程(13)的奇點.接下來我們證明K=K1僅為方程(13)的正則奇點而非本性奇點.

由于tanhK在點K=K1處是解析的，可將tanhK在K=K1處進行Taylor 展開，于是有

因此可得

注意到

不難證明

這證明了K=K1僅為方程(13)的正則奇點而非本性奇點.于是根據(jù)Frobenius 級數(shù)理論[41]，方程(13)的通解可以在K=K1展開成如下Frobenius 級數(shù)：

其中a0(c)=1，而常數(shù)c將由后面出現(xiàn)的指標(biāo)方程決定.將以上通解代入方程(13)，有

其中(K)=(K-K1)S(K)(K)=(K-K1)2T(K).在方程(23)中令K=K1，并利用式(20)和(21)，得到指標(biāo)方程c(c-1)+c/2=0，因此解得c1=0和c2=1/2.通過合并式(23)中K-K1的同次冪項，可以遞推確定系數(shù)am(c)：

其中

由于c1和c2相差一個非整數(shù)，所以方程(13)的兩個特解可分別表示為如下的Frobenius 級數(shù)：

再結(jié)合式(12)，區(qū)間Ij,1上自由液面位移函數(shù)η(x)可表示為以上兩個特解的線性組合：

其中A(1j)和A(2j)為待定系數(shù).

接下來考慮第j個溝槽的后壁區(qū)間Ij,2=[xj,2,xj,3]，水深函數(shù)仍如式(1)所示，但其反函數(shù)與前壁區(qū)間上的反函數(shù)(8)不同，這里變?yōu)?/p>

同樣運用Liu 和Zhou[39]與Liu 和Xie[40]創(chuàng)立的求解修正緩坡方程的變量替換技術(shù)，也在修正緩坡方程(2)中引入變量替換(9)，則由線性波的色散關(guān)系(10)，式(29)可改寫為

注意到水深函數(shù)在Ij,2上也嚴(yán)格單調(diào)遞減，易知由方程(9)和(29)定義的變量替換x→K亦為1-1 映射，它將Ij,2映射到區(qū)間[K0,K1].位移函數(shù)η(x)在變量替換后變?yōu)镵的函數(shù)，記作ζ(K)，即

將式(29)代入文獻[39]中的顯式方程(19)，則第j個溝槽后壁區(qū)間Ij,2上關(guān)于舊變量x的隱式修正緩坡方程(2)被成功轉(zhuǎn)化為區(qū)間[K0,K1]上關(guān)于新變量K的顯式常微分方程：

因為方程(32)與方程(13)在形式上完全一致，故方程(32)的通解也可表為兩個特解ξ1(K)和ξ2(K)的線性組合，所以第j個溝槽后壁區(qū)域Ij,2上的自由液面位移函數(shù)η(x)可表示為

其中A(3j)和A(4j)為待定系數(shù).

為了確定系數(shù)AR，AT以及A(Tj)，j=1,2,···,N，t=1,2,3,4，根據(jù)文獻[42]，自由液面位移函數(shù)η(x)需滿足如下質(zhì)量守恒耦合條件：

具體說來，當(dāng)x=x1,1=0時，條件(34)和(35)變?yōu)?/p>

其中

當(dāng)x=xj,1，j=2,3,···,N時，條件(34)和(35)變?yōu)?/p>

當(dāng)x=xj,2，j=1,2,3,···,N時，條件(34)和(35)變?yōu)?/p>

其中

當(dāng)x=xj,3，j=1,2,3,···,N-1時，條件(34)和(35)變?yōu)?/p>

其中

當(dāng)x=xN,3時，條件(34)和(35)變?yōu)?/p>

最后，結(jié)合式(36)、(40)、(41)、(43)、(45)，可得

因此，由N個周期排列的全等拋物形溝槽引起線性水波反射的反射系數(shù)的解析解為

當(dāng)n=1時，式(46)可簡化為

其中，W(ξ1,ξ2)=ξ1(K0)ξ2′(K0)-ξ2(K0)ξ1′(K0)，|Π0|為復(fù)數(shù) Π0的模，Re(Π0)為復(fù)數(shù) Π0的實部.因此有

2 級數(shù)解收斂的條件以及判斷收斂與否的算例演示

雖然式(28)和(33)給出的是修正緩坡方程的精確解析解，但它們并非封閉解，而只是Frobenius 級數(shù)解.實際計算時，對一列具有特定寬度w、水深h0和h1和間距d的拋物形溝槽，外海入射波周期或頻率的變化，將引起兩個參數(shù)K1和K0的變化.而級數(shù)(26)和(27)并不是在任何情況下都收斂，它們的收斂是有條件的.本節(jié)首先推導(dǎo)出它們收斂的條件，然后通過兩個實際算例來具體演示如何判斷所構(gòu)造的Frobenius 級數(shù)解是否收斂.

首先考察式(26)和(27)中級數(shù)解ξ1(K)和 ξ2(K)的收斂性.注意到方程(13)的所有奇異點如下：

1) sinh(2K)=0的所有復(fù)根，即它們?nèi)课挥谔撦S上，參見圖2 中的紅色圓點.

2) 2K+sinh(2K)=0的所有復(fù)根，參見圖2 中的藍色圓點，其分析可參見文獻[39].

3)KtanhK=K1tanhK1的復(fù)根，它包含兩個實根K=±K1(參見圖2 中的兩個綠色圓點)以及依賴于水深h1和角頻率 ω的一些虛根.因為這些虛根隨h1和 ω變化而變化，所以它們在虛軸上的位置并不固定，因此在圖2中并未標(biāo)出，但我們只需記住它們在虛軸上就可以了.

圖2 方程(13)除K tanh K=K1 tanh K1的虛根外的所有奇異點Fig.2 Singularities of eq.(13) except for the imaginary roots ofK tanh K=K1 tanh K1

根據(jù)Frobenius 級數(shù)理論，兩個級數(shù)ξ1(K)和 ξ2(K)在整個區(qū)間[K0,K1]上收斂的充要條件是，以級數(shù)展開點K1為圓心，且以K1-K0為半徑的復(fù)圓盤|K-K1|≤K1-K0不能包含方程(13)的除K1外的任何奇異點.注意到K0＞0，位于虛軸上的所有奇異點不可能在復(fù)圓盤之內(nèi).而其他不在虛軸上除K1外的所有奇異點中，距離K1最近的奇異點為關(guān)于x軸對稱的兩個奇異點K=1.125 36±2.106 20i.只要這兩個奇異點不落入復(fù)圓盤之內(nèi)，則所有奇異點都不會落入該復(fù)圓盤內(nèi).因此式(26)和(27)中級數(shù)解ξ1(K)和 ξ2(K)收斂的充要條件是如下不等式成立：

上述不等式在幾何上為一個以三條曲線K0=K1-(K1-1.125 36)2+4.436 078 44，K0=K1和K0=0為邊界圍成的右邊沒有邊界的開放區(qū)域，如圖3 所示陰影部分，它即為級數(shù)解ξ1(K)和 ξ2(K)的收斂區(qū)域.當(dāng)且僅當(dāng)點(K1,K0)落在該陰影區(qū)域內(nèi)，ξ1(K)和 ξ2(K)收斂.由圖3 可以看出，當(dāng)K1＜2.6，即溝槽不太深時，對任何K0＜K1，兩個級數(shù)解(26)和(27)都無條件收斂.而當(dāng)2.6 ＜K1＜10，即溝槽適當(dāng)深時，只要K0＞1.009，即海床的整體水深h0不是太淺，則兩個級數(shù)解(26)和(27)也都收斂.

圖3 級數(shù)ξ1(K)和 ξ2(K)的收斂區(qū)域(即陰影區(qū)域)以及兩個例子的收斂性分析Fig.3 Convergent regions of ξ1(K) and ξ2(K) and the convergence analysis in 2 examples

接下來，我們將通過兩個實際的算例來演示如何判斷兩個級數(shù)解(26)和(27)的收斂性.

正如前面所指出的，當(dāng)且僅當(dāng)(K1,K0)落在如圖3 所示的陰影區(qū)域內(nèi)，式(26)和(27)中的級數(shù)解ξ1(K)和ξ2(K)才是收斂的.

進一步地，為了演示如何準(zhǔn)確判斷ξ1(K)和 ξ2(K)收斂，第一個例子我們考慮一組拋物形溝槽，其整體水深h0=0.7m，溝槽最深處水深h1=1.2 m.因為級數(shù)解的收斂性條件(50)與溝槽個數(shù)N、溝槽寬度w和溝槽間距d無關(guān)，因此不用設(shè)定N，w和d的具體數(shù)據(jù).考慮所有滿足條件0.03≤K0≤5.0的入射波，此時入射波周期T的變化范圍為0.75 s≤T≤55.96 s.由線性波色散關(guān)系(10)有

當(dāng)K0由0.03 變化到5.0 時，通過對以上方程數(shù)值求解K1，便可以得出(K1,K0)的移動軌跡如下：

上述(K1,K0)的移動軌跡體現(xiàn)在幾何上為一條曲線線段，在圖3 中以紅色線段標(biāo)注.因為整條線段全部落在陰影區(qū)域中，據(jù)此可以斷定式(26)和(27)中的Frobenius 級數(shù)ξ1(K)和 ξ2(K)均收斂.

第二個例子我們考慮一組拋物形溝槽，整體水深h0=0.2 m，溝槽最深處水深為h1=4.0 m.考慮所有滿足條件0.10≤K0≤0.75的入射波，此時入射波周期T的變化范圍為1.3 s≤T≤8.99 s.類似于第一個算例，可得到K1隨K0變化而變化的如下關(guān)系式：

當(dāng)K0由0.10 變化到0.75 時，對以上方程數(shù)值求解K1，亦可得出(K1,K0)的移動軌跡如下：

上述(K1,K0)的移動軌跡為圖3 中的綠色線段.因該線段沒有完全落在陰影區(qū)域，故ξ1(K)和 ξ2(K)對周期T在1.3 s≤T≤8.99 s中的所有入射波并不能都保證其收斂性.若將K0限制在0.1≤K0≤0.45，也即將周期T限制在2.06 s≤T≤8.99 s，則ξ1(K)和 ξ2(K)都將收斂.

3 模型驗證：與邊界元數(shù)值解的比較

考慮等距排列的拋物形溝槽，其參數(shù)為N=4，h0=1.0 m，w=4 m，d=7 m，然后讓h1變化，分別考慮h1為1 .1 m，1.2m，1.3 m 和 1.4 m 的4 種工況，此時溝槽的無量綱深度H=b/h0分別為0.1，0.2，0.3和 0.4.首先，針對0.002 24≤K0≤2.244，也即0.005 ＜2d/L≤5.0，先采用本文解析模型對這4 種工況進行計算.接下來，為了檢驗本文解析解的正確性，我們采用邊界元數(shù)值模型對0.005≤K0≤2.244，即0.011 14 ＜2d/L≤5.0，也進行了計算，其邊界元模型由滕斌和侯志瑩[43]建立，他們將流域分割為距離溝槽陣列較遠的上游流域、下游流域和溝槽所在的中間流域，在上、下游流域?qū)λ俣葎葑鎏卣髡归_，中間流域則采用邊界元方法離散求解Laplace 方程，最后將三個區(qū)域解聯(lián)立求解.不過，文獻[43]中僅僅計算了地形相對簡單的單個對稱和非對稱二維沙壩和溝槽對波浪的反射，結(jié)果表明在緩坡條件下，邊界元數(shù)值解與修正緩坡方程的解析解[44]吻合得很好.其實，該邊界元模型也可用于計算周期排列的有限系列拋物形溝槽引起的水波Bragg 共振反射，不過地形更復(fù)雜，因而使得處理起來也更復(fù)雜.為了繼續(xù)驗證本文的解析解，我們這里也借用該模型同時計算本節(jié)考慮的4 種工況.

在采用邊界元模型計算時，計算區(qū)域取為[-4(d-w)，Nd+w+4(d-w)]× [-h(x)，0]，在計算區(qū)域的邊界上總共用到1 600 個二次元單元，其中左右邊界均用到50 個二次元單元，而自由面邊界和海底地形邊界均用到750 個二次元單元，其中每個溝槽上50 個單元，兩個相鄰間的溝槽之間50 個單元，溝槽區(qū)域的前后各用到200 個單元.而采用本文解析模型計算時，依據(jù)本文前述級數(shù)解收斂的判斷方法可以很容易判斷出，本文解析模型的級數(shù)解ξ1(K)和 ξ2(K)對于上述4 種工況在整個區(qū)間0.005≤2d/L≤5.0中均收斂.兩種不同模型計算得到的結(jié)果參見圖4.容易看出，對4 組不同深度的溝槽，采用兩種不同模型計算得到的結(jié)果幾乎完全一致，因為4 種情形下最大的海底地形斜率分別為0.1，0.2，0.3 和0.4，都在修正緩坡方程的適用范圍內(nèi)(即海底地形最大斜率不超過1∶1).兩個解相互之間很高的吻合程度佐證了本文解析模型的正確性.

圖4 采用本文解析模型和基于Laplace 方程的邊界元模型[43]分別計算得到的兩個解的比較Fig.4 Comparison between the present solution and the BEM solution[43] to Laplace’s equation

需要指出的是，當(dāng)2d/L很小時，圖4(a)、(b)中邊界元解可見翹尾現(xiàn)象，即當(dāng)2d/L→0時反射系數(shù)并不為0，原因是此時相應(yīng)的入射波波長很長，比如2d/L=0.011 14對應(yīng)的入射波長為L≈1 256.73 m，因此在數(shù)值計算中，若想完整地抓住入射波信息，理論上要求計算區(qū)域在溝槽陣列的前后至少應(yīng)留有一個波長的長度.但是，對于上面考慮的海底系列溝槽的4 個工況，假若溝槽陣列前后真留一個波長的長度，則計算區(qū)域應(yīng)取為[-1 269.23，1 269.23]，而溝槽陣列所在的區(qū)域僅為[-12.5，12.5]，導(dǎo)致溝槽陣列在整個計算區(qū)域中的相對分辨率不夠.因此當(dāng)2d/L很小時，實際計算區(qū)域在溝槽陣列的前后其實又做不到真的留一個波長的長度.圖4(a)～(d)中目前呈現(xiàn)的4 個邊界元解是在溝槽陣列前后各留有4(d-w)=12 m 進行計算的結(jié)果.可見當(dāng)2d/L很小時，邊界元方法中的翹尾現(xiàn)象無法徹底消除.這可能也是很多數(shù)值模型都將計算區(qū)域的左端點取得離0 有一定距離的原因，例如文獻[15]中的圖5、6，文獻[21]中的圖3～7，文獻[45]中的圖2，文獻[46]中的圖1～9.但本文為了邊界元結(jié)果的全域呈現(xiàn)和真實呈現(xiàn)，在采用文獻[43]的邊界元模型計算時仍令2d/L的最小取值為0.011 14，即K0=0.005.

另一方面，本文解析解并不存在任何翹尾的問題.事實上，當(dāng)2d/L=0.005時，本文解析解計算出4 個工況下的反射系數(shù)值分別為0.001 11，0.002 07，0.002 90 和0.003 64.此時入射波波長為2 800 m，4 個寬度僅為4 m，高分別為0.1 m，0.2 m，0.3 m 和0.4 m 的溝槽對于波長如此長的入射波而言，其影響幾乎可以忽略不計，因此從物理上推斷，其反射系數(shù)顯然應(yīng)接近于0，這個推斷與本文解析解計算的結(jié)果完全符合.解析解與邊界元解在這一點上的差異也反映出尋找解析解的重要性，除了計算機位數(shù)引起的舍入誤差(也即機器誤差)外，解析解不存在任何截斷誤差.因此解析解和實驗解都被作為檢驗和驗證數(shù)值解的標(biāo)準(zhǔn)，但解析解的獲取成本卻遠遠低于實驗解的獲取成本.

4 本文解析解與線性攝動解析解的比較

基于經(jīng)典的一階攝動法，Miles[47]建立了連續(xù)微幅波動海底地形或障礙物引起表面波反射的反射系數(shù)的通用解析公式如下：

其中δ(x)為海底地形或障礙物的形函數(shù)，其高度從海底的平均深度h0算起.將該通用解析公式具體應(yīng)用于本文的拋物形系列溝槽，注意到相應(yīng)的形函數(shù)為

可得反射系數(shù)的經(jīng)典攝動解析解如下：

值得指出的是，Miles 的通用公式(55)是基于經(jīng)典的一階攝動法導(dǎo)出的，其適用條件僅限于攝動法的小參數(shù)假設(shè)，即溝槽的無量綱深度H必須很小.另外，從式(57)的第二式可知，當(dāng)沙壩間距d為入射波半波長L/2的整數(shù)倍時，若N足夠大，則反射系數(shù)RMiles有可能無界，說明Miles 的通用公式(55)在N很大時將會崩潰失效.這與Mei[5]指出當(dāng)沙紋個數(shù)足夠大時，Davies[2]關(guān)于正弦沙紋地形Bragg 共振反射的攝動解析解將崩潰失效是完全一致的.

考慮一列共8 個拋物形溝槽，先固定如下參數(shù)：h0=1 m，w=4 m，d=7 m，然后讓h1變化，分別取h1為1.05 m，1.15m，1.25 m 和1.35 m，此時溝槽的無量綱深度H分別為0.05，0.15，0.25和0.35.分別采用本文基于修正緩坡方程的級數(shù)形式解析解模型和Miles 基于經(jīng)典一階攝動法的封閉形式解析解模型對此進行計算，計算結(jié)果參見圖5.如圖5(a)所示，當(dāng)H=0.05時，采用兩種方法計算得到的結(jié)果幾乎完全一致，這就互相印證了兩個結(jié)果的正確.但隨著H的值逐步增大，即當(dāng)H=0.15和0.25 時，兩個解的差別也隨之逐步顯現(xiàn)，參見圖5(b)、(c)，原因在于Miles 基于經(jīng)典一階攝動法的通用公式(55)已經(jīng)不再適用，因而產(chǎn)生了較大誤差.最后，當(dāng)H=0.35時，基于一階攝動法的解析解模型完全崩潰，因為采用它所計算出來的反射系數(shù)竟然超過了1，達到了1.115 3，這顯然是錯誤和荒謬的.

圖5 本文解析解與應(yīng)用Miles 的通用公式(55)的比較Fig.5 Comparison between the present solution and the solution based on Miles’ formula (55)

作為對比，本文基于修正緩坡方程建立的解析模型并不受小參數(shù)的限制，唯一要求的是海底地形變化相對較緩.正如文獻[44]所顯示的，在計算單個梯形沙壩或溝槽引起的波浪反射時，只要海底地形斜率不超過1∶1，則基于修正緩坡方程的解析模型與文獻[48]的基于Laplace 方程的特征函數(shù)展開模型(EFEM)的解吻合得很好，參見文獻[44]中的圖3(b)、5(b)和6(b).因此本文建立的基于修正緩坡方程的解析解與Miles 的攝動解析解相比，其適用范圍大，既不受溝槽個數(shù)的限制，也適用斜率不超過1∶1 的溝槽.

5 溝槽個數(shù)對Bragg 共振反射峰值、帶寬和相位的影響

為研究溝槽個數(shù)對反射系數(shù)特別是共振反射系數(shù)的影響，考察具如下參數(shù)的系列溝槽：h0=0.7 m，h1=1.2 m，w=2.0m，d=3.0 m，并令溝槽個數(shù)N從1 增加到20，計算范圍為0.05≤2d/L≤3.0，等價于2 m ≤L≤ 120 m.由第2 節(jié)討論知，本文解析模型中的級數(shù)解ξ1(K)和 ξ2(K)對所有20 種情形在整個范圍0.05≤2d/L≤3.0中均收斂.相應(yīng)計算結(jié)果參見圖6.

可以看出，當(dāng)溝槽個數(shù)N增加時，Bragg 共振反射的所有主振峰值(圖6 中以圓點表示)均相應(yīng)地逐步增加，很快就達到了1.而相應(yīng)共振帶的帶寬則隨溝槽個數(shù)增加而逐步變窄，依據(jù)文獻[13, 49]，Bragg 共振反射的共振帶以無限周期列拋物形溝槽相應(yīng)的Bloch 波帶隙(band gap)為極限，這與有限周期排列人工沙壩導(dǎo)致的現(xiàn)象一致，參見文獻[38]的圖3 和圖11.

圖6 拋物形溝槽個數(shù)對Bragg 共振反射及共振帶寬的影響Fig.6 Influences of the number of trenches on the Bragg resonance and the resonance bandwidth

現(xiàn)在來看Bragg 共振主振發(fā)生的相位，它們均向右發(fā)生了偏移，即發(fā)生在2d/L＞1的位置.在此之前，不少的計算結(jié)果，包括數(shù)值解[28,49]、實驗解[28]和解析解[11,24-26]均表明，由海底周期排列的正弦沙紋地形或人工沙壩導(dǎo)致的海洋表面波的Bragg 共振反射[13]，并不遵循傳統(tǒng)的Bragg 原理，精確地出現(xiàn)在2d/L=1，而是出現(xiàn)在2d/L＜1的位置，簡稱相位下移(phase downshift)，參見文獻[8-9, 11].最近，Liu 等[14]建立的基于修正緩坡方程的解析模擬表明，由正弦沙紋地形引起的Bragg 共振，若沙紋個數(shù)固定，令沙紋高度增加時，相位下移現(xiàn)象會變得更加明顯，參見文獻[14]中的圖11(c)～(k).Liu 等針對三種特定沙紋地形給出的135 個工況的計算結(jié)果進一步揭示出，當(dāng)正弦沙紋個數(shù)大于8 時，Bragg 主振峰的相位將不再隨沙紋個數(shù)的改變而改變，而是固定不動，參見文獻[14]中的表2.當(dāng)正弦沙紋個數(shù)大于12 時，Bragg 高階共振峰的相位也不再隨沙紋個數(shù)的改變而改變，也是固定不動，參見文獻[14]中的表3.關(guān)于人工沙壩系列引起的Bragg 共振相位下移的研究參見文獻[24-26].類似地，幾乎所有這些早期研究的零星算例(數(shù)值的或?qū)嶒灥?都表明，當(dāng)人工沙壩個數(shù)改變時，相位下移的程度也隨之改變.而Liu 等[30]針對一類特定梯形沙壩的30 個工況的計算結(jié)果則進一步揭示出，前述結(jié)論僅對沙壩個數(shù)很少的情形成立，一旦沙壩個數(shù)相對比較大時，共振峰值的相位下移將不再隨沙壩個數(shù)的改變而改變，而是固定不動，參見文獻[30]中的表1.

回到上面所說的Bragg 共振相位向右邊高頻方向移動的現(xiàn)象.該現(xiàn)象最近才剛剛由Guo 等[35]第一次觀測到.他們研究的是周期排列的旋輪線形系列溝槽所激發(fā)的水波Bragg 共振反射，與以前發(fā)現(xiàn)的相位下移現(xiàn)象相對照，此現(xiàn)象被Guo 等[35]命名為“相位上移”(phase upshift).本文中，相位上移現(xiàn)象再次在周期排列的拋物形系列溝槽上得到確認(rèn)，說明了針對周期排列的溝槽地形，它極可能是一個普遍現(xiàn)象.換言之，只要是周期排列的系列溝槽，無論溝槽形狀如何(如常見的矩形、梯形、三角形、半余弦、半圓形等)，相應(yīng)的共振相位就會上移.我們預(yù)測，這個最近才剛剛發(fā)現(xiàn)的新現(xiàn)象后續(xù)有望吸引學(xué)者們開展進一步的深入研究.

對于由正弦沙紋和人工沙壩引起的Bragg 共振反射出現(xiàn)的相位下移現(xiàn)象，雖然背后隱藏的定量化數(shù)學(xué)機理尚不清楚，但關(guān)于其定性化的物理機制的解釋已經(jīng)不少.Guazzelli 等[9]將此現(xiàn)象的發(fā)生歸于瞬逝波模態(tài).Liu 和Yue[11]將此現(xiàn)象的發(fā)生歸于表面波的非線性性(波陡)，海底地形變化涉及的非線性效應(yīng)(坡陡)，以及水深條件等因素不能嚴(yán)格滿足線性理論的假設(shè).Chang 和Liou[24]將此現(xiàn)象的發(fā)生歸于沙壩對反射波相速度的減緩.最近，Liang 等[18]基于Mathieu 不穩(wěn)定理論給出了波浪在正弦沙紋地形上傳播時Bragg 共振相位的理論預(yù)測公式.

本文研究的是線性波在系列溝槽上的傳播，并非系列沙壩上的傳播，且沒有考慮瞬逝波模態(tài)，因此Bragg共振相位上移現(xiàn)象可解釋為周期排列溝槽的存在使得水波傳播速度加快，因為波浪傳播相速度近似與水深的平方根成正比.事實上，對近岸淺水波，相速度可近似簡化為可見水越深水波傳播速度越快.下面以兩個溝槽為例來解釋它們所激發(fā)的Bragg 共振反射相位上移現(xiàn)象.

首先，假定溝槽非常淺非常淺，根據(jù)傳統(tǒng)Bragg 原理，Bragg 共振反射發(fā)生的條件為：2d=nL，其中n為正整數(shù).其物理解釋為，入射波沿x軸正向傳播遇到第1 個溝槽時，將產(chǎn)生沿x軸負(fù)向傳播的反射波和繼續(xù)正向傳播的透射波，為敘述方便，分別稱它們?yōu)榈? 反射波和第1 透射波.第1 透射波繼續(xù)正向傳播距離d后，在遇到第2 個溝槽時又將再次產(chǎn)生負(fù)向傳播的反射波和繼續(xù)正向傳播的透射波，分別稱它們?yōu)榈? 反射波和第2 透射波.而第2 反射波負(fù)向回傳距離d后，再經(jīng)過第一溝槽時又將產(chǎn)生負(fù)向傳播的透射波和正向傳播的反射波，為描述的統(tǒng)一起見，稱其中負(fù)向傳播的透射波為第3 反射波，而正向傳播的反射波為第3 透射波.所謂Bragg 共振反射，指的是第1 反射波與第3 反射波正好同相位疊加.由于此時假定兩個溝槽都非常淺，因此無論是入射波、反射波，還是透射波，其傳播的相速度都可以看做近似相同，記作C.因此，要第1 反射波與第3 反射波同相位疊加產(chǎn)生Bragg 共振反射，則第3 反射波自離開第一溝槽到再次回到第一溝槽期間走過距離 2d所用的時間，必須恰好等于第1 反射波走過其波長L的整數(shù)倍n所花的時間，即2d/C=nL/C，也即2d/L=n.

現(xiàn)在假定考慮的兩個拋物形溝槽具有一定深度，雖然第3 反射波自離開第一溝槽繼續(xù)正向傳播到被第二溝槽反射后再次回傳并穿過第一溝槽所走過的距離仍然為 2d，但是由于溝槽具有一定深度，第3 反射波在通過這個距離 2d時的相速度比原來的要快，假定其平均相速度為C′，則C′＞C.因此，要第1 反射波與第3 反射波同相位疊加產(chǎn)生Bragg 共振反射，則第3 反射波自離開第一溝槽到再次回傳后穿過第一溝槽所走過距離 2d所用的時間，必須恰好等于第1 反射波走過其波長L的整數(shù)倍n所 花的時間，即2d/C′=nL/C，也即2d/L=C′n/C＞n.顯然，此時Bragg 共振反射出現(xiàn)了相位上移現(xiàn)象.

6 溝槽深度對Bragg 共振反射峰值、帶寬和相位的影響

為討論拋物形溝槽深度b對Bragg 共振反射的影響，對如下固定的溝槽參數(shù)h0=0.5 m，w=2.0 m，d=3.0 m，N= 4，6，8，10，令溝槽的無量綱深度H取如下5 個值：0.2，0.3，0.4，0.5，0.6，也即b= 0.10 m，0.15 m，0.20 m，0.25 m，0.30 m，然后采用本文建立的基于修正緩坡方程的解析模型計算了當(dāng)0.05≤2d/L≤3.0，亦即2 m≤L≤120 m 時的反射系數(shù)，其結(jié)果參見圖7.

從圖7 可以看出，當(dāng)溝槽個數(shù)N、寬度w及相鄰溝槽的間距d都固定的情況下，Bragg 主振峰值隨溝槽深度的增加而增加.其次，當(dāng)溝槽深度增加時，共振帶寬也隨之增加但增加的幅度不大.最后還可看出，在溝槽個數(shù)N、寬度w以及間距d固定的情況下，共振相位隨溝槽深度的增加向高頻移動的幅度也增加，其原因與前面解釋的一致，即溝槽越深，第3 反射波自離開第一溝槽繼續(xù)正向傳播到被第二溝槽反射后再次回傳并穿過第一溝槽的平均相速度更快，從而導(dǎo)致共振相位的上移幅度更大.

7 溝槽寬度對Bragg 共振反射峰值、帶寬和相位的影響

最后，我們討論溝槽寬度對Bragg 共振反射峰值、帶寬和相位的影響.對如下固定的溝槽參數(shù)：h0=0.7 m，b=0.5m (即h1=1.2 m)，d=12 m，N= 4，然后令寬度w分別取1 m，2 m，···，12 m.采用本文建立的基于修正緩坡方程的解析模型對0.05≤2d/L≤3.0，亦即8 m≤L≤ 480 m 時的反射系數(shù)進行計算，其結(jié)果參見圖8(a)～(c).

圖8 溝槽寬度對Bragg 共振的影響Fig.8 Influences of the trench width on the Bragg resonance

先看溝槽寬度對Bragg 共振反射峰值的影響.將圖8(a)～(c)中標(biāo)記出的主振反射峰值以及相對應(yīng)的相位整合至圖8(d)中，可以清晰地看出當(dāng)拋物線形溝槽的寬度w由1 m 增加到8 m 時，Bragg 共振反射的峰值由0.159 18單調(diào)增加至0.729 48，但當(dāng)溝槽寬度w由8 m 繼續(xù)增加時，Bragg 共振反射的峰值則反而開始逐步降低，一直到w=12 m 時，Bragg 共振反射的峰值降為0.594 97.這就意味著在這個特殊的工況下，Bragg 共振反射的峰值在w=8m 的附近達到最大值.其次還可以看出，拋物線形溝槽的寬度w變化時，Bragg 共振帶寬幾乎不受影響.最后，在拋物線形溝槽個數(shù)N、溝槽深度b以及相鄰溝槽之間間距d固定的情況下，Bragg 共振反射的主振相位隨溝槽寬度的增加向高頻移動的幅度也增加，其原因也與前面解釋的一致，即溝槽越寬，水變深的范圍更大，第3 反射波自離開第一溝槽繼續(xù)正向傳播到被第二溝槽反射后再次回傳并穿過第一溝槽的平均相速度更快，從而也導(dǎo)致共振相位的上移幅度更大.

8 結(jié) 論

針對周期排列的拋物線形有限系列溝槽對線性表面波的反射問題，本文構(gòu)造了修正緩坡方程的Frobenius 級數(shù)解.在此基礎(chǔ)上，結(jié)合矩陣乘法，給出了反射系數(shù)的解析公式.為了檢驗解析解的正確性，我們將其與邊界元解進行了對比，發(fā)現(xiàn)二者吻合很好，但邊界元解對波長非常長的入射波有翹尾現(xiàn)象，意味著一定誤差的存在，而本文解析解沒有翹尾現(xiàn)象.與Miles 的一階攝動解的比較反映了本文解析解適用性范圍更廣，可以適用于深度較大的拋物溝槽(只要溝槽最大斜率不超過1∶1)，而攝動解在溝槽地形斜率很大時可能崩潰失效.

基于本文建立的反射系數(shù)的解析公式，我們分別討論了拋物形溝槽個數(shù)、深度和寬度對Bragg 共振反射峰值、帶寬和相位的影響.結(jié)果表明，當(dāng)固定溝槽寬度與深度而令溝槽個數(shù)增加時，Bragg 共振峰值將隨之相應(yīng)地增大，直到最終達到1，即全反射，而Bragg 共振帶寬則逐步變窄，終止于一個固定的寬度而不再進一步變窄.當(dāng)固定溝槽個數(shù)和溝槽寬度，而令溝槽深度增加時，Bragg 共振峰值將隨之相應(yīng)地逐步增大，共振帶寬也隨之稍微增加但增加的幅度不大.當(dāng)固定溝槽個數(shù)和溝槽深度，而令溝槽寬度增加時，Bragg 共振峰值則先增后減，說明存在一個特別的溝槽寬度，使得Bragg 共振反射峰值達到最大.當(dāng)固定溝槽個數(shù)和溝槽深度，而令溝槽寬度增加時，Bragg 共振帶寬基本不變.最重要的是，本文確認(rèn)了周期排列的有限個拋物形溝槽也將導(dǎo)致Bragg 共振的相位上移.該現(xiàn)象的再次確認(rèn)說明了，相位上移對于有限個周期排列的任意形狀的溝槽地形極大可能是一個普遍現(xiàn)象，原因在于任意形狀溝槽的存在都使得水深變深了.另外，計算結(jié)果顯示，當(dāng)增加溝槽深度或?qū)挾葧r，相位上移的幅度也隨之增大，因為水深變深的程度隨之增大.

致謝本文作者衷心感謝浙江海洋大學(xué)科研啟動基金(Q1607)對本文的資助和大連理工大學(xué)滕斌教授提供的有關(guān)圖4 中邊界元數(shù)值結(jié)果的計算程序.