琚海軍，劉忠波，劉 勇，房克照

(1.中國海洋大學(xué) 工程學(xué)院，青島 266100；2.大連海事大學(xué) 交通運(yùn)輸工程學(xué)院，大連 116026；3.大連理工大學(xué) 海岸和近海工程國家重點(diǎn)實驗室，大連 116023)

為了保護(hù)海岸，國內(nèi)外修建了不少帶有多孔介質(zhì)的拋石防波堤等海岸工程。這些建筑物本身帶有很大孔隙，它們一方面可以將遠(yuǎn)海傳來的波浪反射回去，另一方面也可通過建筑物自身內(nèi)部的摩擦阻力消耗掉一部分波浪，對波浪運(yùn)動產(chǎn)生明顯的衰減影響，從而起到了保護(hù)海岸的效果。考慮到近岸海域的海床上覆蓋著大片透水程度不同的顆粒介質(zhì)，這些實際孔隙介質(zhì)的滲透性也直接影響到岸灘上的波浪水流等水動力特性，所以滲透性對波浪運(yùn)動產(chǎn)生的影響也是不能忽略的。

關(guān)于孔隙海床的波浪水動力特性研究，多集中于孔隙海床上波高衰減問題上，并且常借助于實驗或理論分析等手段進(jìn)行展開。過去一些研究者對滲透海床上的波浪水動力問題進(jìn)行了初步探究和分析，但是這些模型不能直接應(yīng)用于研究較大范圍內(nèi)可滲海床上的波浪信息，而模擬近岸海域附近的數(shù)學(xué)模型主要有基于緩坡假定的緩坡方程和基于小振幅假定的Boussinesq型水波方程。緩坡方程的模型則多以線性波浪為研究對象且僅能適應(yīng)于考慮弱非線性的情況，無法用于非線性較強(qiáng)的波浪與波浪、波浪與水流等相互作用問題的研究。而高階Boussinesq型水波方程在模擬強(qiáng)非線性波浪、波浪與波浪相互作用等相關(guān)問題上則顯示出較好的適用性。

但是國內(nèi)外研究者在推導(dǎo)這些方程的時候，為了簡化問題，通常假設(shè)海底不透水[1-7]，僅有少數(shù)研究者進(jìn)行可滲海床上這方面的研究[8-13]。而后一類模型也只能考慮上部水體為自由水體，下部為可滲多孔介質(zhì)的波浪運(yùn)動情況，這使得它們很難直接應(yīng)用于模擬波浪在多孔介質(zhì)中傳播變形。Hsiao等[12]通過歐拉方程給出了適應(yīng)于多孔介質(zhì)中波浪傳播變形的高階Boussinesq型水波方程，模型是由可滲透多孔介質(zhì)某一深度處速度和波面升高表達(dá)的。劉忠波等推導(dǎo)了兩層滲透介質(zhì)上的Boussinesq型水波方程[13]，我們通過對其一組方程進(jìn)行了簡化，我們導(dǎo)出了一組高階Boussinesq型水波方程，對其進(jìn)行了理論研究，并在此基礎(chǔ)上進(jìn)行了數(shù)值研究工作。

1 高階Boussinesq型水波方程的理論推導(dǎo)

波浪在多孔介質(zhì)中的傳播見圖1，η(x,y,t)表示波浪在多孔介質(zhì)中的波面升高(以靜水面為起算面)，多孔介質(zhì)的厚度為h(x,y)，ψ(x,y,z,t)表示速度勢。

圖1 波浪在多孔介質(zhì)中傳播的示意圖Fig.1 Sketch of wave propagation in porous media

在無因次化中常用的三個特征量為：a、h0和l0，為了確認(rèn)方程中不同項的相對重要性，引入以下無因次

(1)

在Boussinesq類水波方程推導(dǎo)過程中，常引入非線性特征參數(shù)和色散性特征參數(shù)

(2)

為了簡便起見，忽略了上標(biāo)，則對應(yīng)的控制方程和邊界條件可寫為

μ2▽2ψ+ψzz=0，-hb

(3)

(4)

μ2(ηt+ε▽ψ·▽η)=ψz，z=εη

(5)

μ2▽ψ·▽h+ψz=0，z=-h

(6)

這里α=α1+α2|U|，其中α1和α2分別代表層流和紊流阻力系數(shù)；cr=1+(1-n)cm/n為慣性系數(shù)，cm為附加質(zhì)量系數(shù)，詳細(xì)可參見文獻(xiàn)[8-9]。

將速度勢展開如下

(7)

分別對速度勢求1次和2次導(dǎo)數(shù)，并將結(jié)果帶入控制方程(3)中可得

(8)

所以速度勢表達(dá)式可以表示為

(9)

將(9)代入(6)中，經(jīng)整理可得

(10)

式中：w0和u0分別表示在靜止水位初的波浪垂向和水平速度。

對方程(4)求導(dǎo)，并將(10)代入，整理可得

(11)

類似地，整理(4)得到的連續(xù)方程如下

ηt+▽·Q=0

(12)

引入水深積分平均速度

(13)

對(13)進(jìn)行整理，可得

(14)

連續(xù)方程可寫為

(15)

將(14)代入(11)中，整理可得

(16)

方程(15)和方程(16)組成了一組Boussinesq方程。

沿水深積分平均速度表達(dá)的方程色散關(guān)系式與經(jīng)典Boussinesq方程的一致，因此該方程的適用水深有限。為了拓展方程的適用水深，在動量方程(16)中引入劉忠波等[7]建議的表達(dá)式

(17)

式(16)中的兩個參數(shù)可以改善方程的色散性，當(dāng)增加(17)后，方程(16)可以寫成

(18)

由(15)和(18)組成的方程，當(dāng)不考慮α和cr時，方程轉(zhuǎn)化為劉忠波等給出的適合非滲透Boussinesq水波方程[7]。同時，當(dāng)文獻(xiàn)[13]不考慮第二層滲透時，也可以轉(zhuǎn)化為本文的研究方程。

2 方程的色散性與衰減率特征分析

在不考慮海底地形變化情況下，忽略非線性和忽略水深變量對x的導(dǎo)數(shù)，對由方程(15)和(18)組成的模型進(jìn)行色散性分析，可分別得到如下色散關(guān)系式

(19)

式中:ω代表頻率；Kb為復(fù)波數(shù)，Kb=Krb+iKib，虛數(shù)i滿足i2=-1,β=β1+β2。

為了對比不同方程的相速度和ki，采用了Hsiao等[12]文獻(xiàn)中給出的解析色散關(guān)系表達(dá)式

ω2(cr+α1i/ω)=gKtanh(Kh)

(20)

這里K=krs+ikis為復(fù)波數(shù)。

對比(19)和(20)，可以測定不同情況下方程的相速度與衰減率，反推參數(shù)的取值。在這里, 從0.1 s/m變化到4 s/m，krh從0.1變化到3，通過比較，以下面兩種判斷式為標(biāo)準(zhǔn)我們給出第一組解和第二組解。

[1-(krs/krb)]2+[1-(kis/kib)]2=min (21)

由于方程中有兩個參數(shù)，而根據(jù)以上對比，只能確定一個參數(shù)，因此，另外一個參數(shù)則只能通過非滲透性情況下的變淺性能獲取，具體可參照文獻(xiàn)[7]，這里直接給出結(jié)果，見表1。對應(yīng)的無因次相速度、衰減率以及復(fù)波數(shù)見圖2和圖3。從衰減率和相速度來看，第一組參數(shù)的適用范圍看略大一些。從復(fù)波數(shù)變化來看，第二組的適用范圍更大。整體來看，二者差異性并不十分明顯，在后文的數(shù)值模擬中采用了第二組參數(shù)。

2-a 無因次相速度2-b 無因次衰減率2-c 無因次復(fù)波數(shù)圖2 方程的無因次相速度(a)、無因次衰減率(b)和無因次復(fù)波數(shù)(c)(第一組)Fig. 2 The phase velocity (a), the dimensionless attenuation rate (b), and the dimensionless complex wave number (c) of the equation. (Group 1)

3-a 無因次相速度3-b 無因次衰減率3-c 無因次復(fù)波數(shù)圖3 方程的無因次相速度(a)、無因次衰減率(b)和無因次復(fù)波數(shù)(c)(第二組)Fig. 3 The phase velocity (a), the dimensionless attenuation rate (b), and the dimensionless complex wave number (c) of the equation. (Group 2)

3 算例

3.1 孤立波穿越多孔直立防波堤

針對上面給出的模型在非交錯網(wǎng)格下進(jìn)行了離散，時間步進(jìn)采用了三階Adams-Bashforth格式進(jìn)行預(yù)報，進(jìn)而采用四階的Adams-Moulton格式進(jìn)行校正[7]。當(dāng)預(yù)報值和校正值滿足一定誤差時，當(dāng)前步結(jié)束，否則重新更新變量，重新利用校正步進(jìn)行校正。

圖4 實驗布置Fig.4 Experimental arrangement

用本文提出的理論模擬Vidal等進(jìn)行的孤立波與多孔直立式防波堤之間的水槽實驗[14]，測量反射系數(shù)和透射系數(shù)。在實驗中，防波堤的的寬度為20 cm或40 cm，實驗水深從25 cm到32 cm不等，防波堤材料的孔隙率為n=0.44，粒徑采用兩種d50=1.43 cm和2.43 cm，實驗中采用的兩個浪高儀分別安置在防波堤前3.9 m和后2 m的地方，如圖4。

在數(shù)值模擬中，模型采用了空間步長△x= 0.01 m，時間步長為△t=0.005 s，圖5給出了本文模型模擬結(jié)果與試驗結(jié)果的比較。在本文計算中，針對d50=1.43 cm，參數(shù)α1和α2取值為18 s-1和234.45 m；針對d50=2.43 cm，參數(shù)α1和α2取值為6.23 s-1和137.97 m，其中α1的值參照Gent給出的系數(shù)確定[15]，α2根據(jù)Lin和Karunarathna給出的系數(shù)確定[16]。由圖5可見，其反射系數(shù)和透射系數(shù)的數(shù)值結(jié)果與試驗結(jié)果相對吻合。

5-a a=20 cm,h=30 cm,d50=1.43 cm5-b a=20 cm,h=30.2 cm,d50=2.43 cm5-c a=40 cm,h=30 cm,d50=1.43 cm5-d a=40 cm,h=31.7 cm,d50=2.43 cm注:符號為Vidal等(1988)試驗數(shù)據(jù);實線:數(shù)值反射系數(shù);虛線:數(shù)值透射系數(shù)圖5 孤立波與多孔直立式防波堤相互作用的反射系數(shù)和透射系數(shù)比較Fig.5 Comparison of reflection and transmission coefficients of solitary wave interaction with porous breakwaters

圖6 實驗布置Fig.6 Experimental arrangement

3.2 滲透地形的數(shù)值實驗

與文獻(xiàn)[17]物理模型實驗類似，僅在8～10 m(堤頂)上面放置出水滲透介質(zhì)。在這里只考慮0.4 m水深，波高為0.02 m，周期T=2.02 s。水槽長23 m、寬0.8 m、高0.8 m，水槽中設(shè)置一梯形障礙物(設(shè)比尺為l:50時，對應(yīng)實際地形中的沙壩)。梯形頂端長2.0 m，頂端的中心位置距離水槽左端為9.0 m。水深為0.4 m，其實驗裝置如圖6。

…… (改進(jìn)的Boussinesq模型) ——— (對照模型)

圖7 改進(jìn)Boussinesq方程數(shù)值波面時間歷程圖對比Fig.7 Comparison of the numerical wave surface of the improved Boussinesq equation

假設(shè)本文方程參數(shù)α2=0，以α1取2.0的計算(以α1=0為對照模型)，T= 2.02 s，附加水體質(zhì)量為0.34，孔隙率為n=0.5，得出對比結(jié)果。圖7對應(yīng)參數(shù)α1取2.0時各個位置的波形圖。通過與對照模型的比較，其結(jié)果驗證了改進(jìn)的Boussinesq型水波方程的有效性，波浪衰減非常明顯，反映出滲透對波浪的影響。

圖8和圖9是α1分別等于0.2，0.5，1.0，2.0的均方波高(統(tǒng)計最后一個周期)和t=40 s時的空間波形比較的結(jié)果。我們也可以明顯得到參數(shù)α1會影響波浪的衰減，并且在同等條件下，隨著α1的系數(shù)值增加，波高衰減增大。

圖8 均方波高(對照組:α1=0)Fig.8 Mean square wave height (control group:α1=0)圖9 t=40 s時的空間波形(對照組:α1=0)Fig.9 Spatial waveform (t= 40 s,control group:α1=0)

4 結(jié)論

在考慮多孔拖曳阻力和慣性阻力的基礎(chǔ)上，推導(dǎo)了一組高階Boussinesq型水波方程。在常水深情況下，分析了模型對應(yīng)的色散關(guān)系式，并與解析解進(jìn)行了比較。建立了Boussinesq型水波方程數(shù)值模型，并模擬了孤立波穿越滲透多孔直立式防波堤和波浪在滲透地形下傳播。通過以上研究，得到結(jié)果如下：

(1)Boussinesq型水波方程中的參數(shù)標(biāo)準(zhǔn)不同時，將給出不同的參數(shù)值；其中以復(fù)波數(shù)比為標(biāo)準(zhǔn)下，得到的參數(shù)與文獻(xiàn)[7]一致。

(2)本文推導(dǎo)的Boussinesq型水波方程可以看成是對文獻(xiàn)[7]的拓展，相當(dāng)于在原方程中多考慮多孔拖曳阻力和慣性阻力。

(3)數(shù)值結(jié)果與實驗結(jié)果吻合良好，說明本文建立的Boussinesq型水波方程是有效的。