滲透介質中波浪傳播的二階Boussinesq水波數值模擬

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

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

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

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

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

1 高階Boussinesq型水波方程的理論推導

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

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

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

(1)

在Boussinesq類水波方程推導過程中，常引入非線性特征參數和色散性特征參數

(2)

為了簡便起見，忽略了上標，則對應的控制方程和邊界條件可寫為

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

(3)

(4)

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

(5)

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

(6)

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

將速度勢展開如下

(7)

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

(8)

所以速度勢表達式可以表示為

(9)

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

(10)

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

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

(11)

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

ηt+▽·Q=0

(12)

引入水深積分平均速度

(13)

對(13)進行整理，可得

(14)

連續方程可寫為

(15)

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

(16)

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

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

(17)

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

(18)

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

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

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

(19)

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

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

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

(20)

這里K=krs+ikis為復波數。

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

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

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

2-a 無因次相速度2-b 無因次衰減率2-c 無因次復波數圖2 方程的無因次相速度(a)、無因次衰減率(b)和無因次復波數(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 無因次復波數圖3 方程的無因次相速度(a)、無因次衰減率(b)和無因次復波數(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 孤立波穿越多孔直立防波堤

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

圖4 實驗布置Fig.4 Experimental arrangement

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

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

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)試驗數據;實線:數值反射系數;虛線:數值透射系數圖5 孤立波與多孔直立式防波堤相互作用的反射系數和透射系數比較Fig.5 Comparison of reflection and transmission coefficients of solitary wave interaction with porous breakwaters

圖6 實驗布置Fig.6 Experimental arrangement

3.2 滲透地形的數值實驗

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

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

圖7 改進Boussinesq方程數值波面時間歷程圖對比Fig.7 Comparison of the numerical wave surface of the improved Boussinesq equation

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

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

圖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 結論

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

(1)Boussinesq型水波方程中的參數標準不同時，將給出不同的參數值；其中以復波數比為標準下，得到的參數與文獻[7]一致。

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

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