劉忠波，孫昭晨，房克照

（大連理工大學(xué) 海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024）

0 引 言

波浪是近岸海域重要的水動(dòng)力因素，準(zhǔn)確掌握波況對(duì)海岸工程建設(shè)有非常重要的現(xiàn)實(shí)意義.研究波浪傳播變形通常有現(xiàn)場(chǎng)觀測(cè)、室內(nèi)物理模型實(shí)驗(yàn)或計(jì)算機(jī)數(shù)值模擬等方法，相較而言，數(shù)值模擬省時(shí)、省費(fèi)用，但其前提是存在較為合適的波浪數(shù)學(xué)模型.近年來(lái)，作為模擬非線性波浪運(yùn)動(dòng)的一類(lèi)高效的數(shù)值模型，Boussinesq類(lèi)水波方程發(fā)展很快［1－3］，但是大多模型是在假定海床不透水的情況下給出的，不能直接用來(lái)模擬波浪在由砂石、砂礫組成的可滲海床上的傳播變形.Cruz等第一個(gè)給出了適合滲透海床上的改進(jìn)型Boussinesq方程［4］，該方程由于僅含弱非線性項(xiàng)，難以描述強(qiáng)非線性波浪運(yùn)動(dòng).Hsiao等、Chen 從歐拉方程出發(fā)，分別給出具有二階色散和二階非線性的高階Boussinesq方程［5－6］，其中Hsiao等也給出了一組以沿深度積分平均速度表達(dá)的Boussinesq 方程［5］，雖然其含有二階非線性項(xiàng)，但其色散性與經(jīng)典Boussinesq方程一致［7］，限制了其在較深水域的應(yīng)用.為了改善該方程的色散性能，劉忠波等增加了二階色散項(xiàng)［8］，并從理論方面分析和討論了不同參數(shù)值下相速度及衰減率的變化.Avgeris等直接在二階非線性Boussinesq方程中耦合可考慮滲透介質(zhì)存在時(shí)的阻力方程，建立了適用于可滲海床情況的Boussinesq類(lèi)水波方程，并進(jìn)行了波浪在可滲潛堤上的數(shù)值研究，數(shù)值結(jié)果表明，這一方式是有效的［9］.而這種方法對(duì)具有四階色散精度的Boussinesq方程（適合非滲透情況）是否有效，目前尚未有相關(guān)的研究結(jié)果發(fā)表.本文針對(duì)這一問(wèn)題進(jìn)行探討，在文獻(xiàn)［1］的基礎(chǔ)上，嘗試引入Cruz等的動(dòng)量方程［4］，聯(lián)立組成適用于可滲海床的Boussinesq方程，從理論上分析方程的相速度以及衰減率，進(jìn)而建立相關(guān)數(shù)值模型，并利用數(shù)值模擬波浪在滲透潛堤和滲透地形上傳播變形，最后通過(guò)相關(guān)實(shí)驗(yàn)結(jié)果檢驗(yàn)該模型的適用程度.

1 基本方程

劉忠波等推導(dǎo)了一組具有四階色散精度的Boussinesq類(lèi)水波方程［1］，該方程表達(dá)如下：

式中：下標(biāo)t表示變量對(duì)時(shí)間的偏導(dǎo)數(shù)為水平二維偏微分算子，η為波面升高，U為沿水深積分的平均速度，h為靜水面以下水深，g為重力加速度.

方程（1）和（2）組成了一組高階Boussinesq水波方程，但該方程是基于海底不可滲這一假設(shè)下給出的，不能應(yīng)用于可滲海床情況.為了考慮滲透海床對(duì)波浪運(yùn)動(dòng)產(chǎn)生的影響，須耦合考慮滲透海床內(nèi)流體運(yùn)動(dòng).具體做法如下：首先保證新建方程的連續(xù)方程是守恒的，因此在連續(xù)方程中引入n·（hsUs）項(xiàng)（Us為滲透介質(zhì)中沿水深積分的平均速度，n為孔隙率）；其次是需引入滲透介質(zhì)中流體運(yùn)動(dòng)方程，本文采用Cruz等的結(jié)果［4］；最后在自由水體中的動(dòng)量方程中需考慮滲透耦合影響，這里也直接采用Cruz等的結(jié)果，也就是引入.最終方程具有以下形式：

式中：hb＝h＋hs，h和hs分別為上層水體的靜水深和滲透介質(zhì)的厚度；β為色散系數(shù)，β＝1／15；cr＝1＋cm（1－n），cm為附加質(zhì)量系數(shù)；α＝a1＋分別為線性阻力系數(shù)和非線性阻力系數(shù).方程（3）～（5）組成了一組適合滲透地形的高階Boussinesq水波方程.

2 理論分析

為了解方程的色散性能，在一維情況下對(duì)方程的色散關(guān)系進(jìn)行理論分析.不考慮水深和滲透介質(zhì)厚度對(duì)空間的導(dǎo)數(shù)，同時(shí)忽略方程中所有非線性項(xiàng)，方程（3）、（4）和（5）可寫(xiě)為

式中：r＝hs／h，將（η，u，us）＝（η0，u0，us0）ei（kx－ωt）代入上述方程中，可以得到

式 中：A11＝－ω；A12＝kh；A13＝nrkh；A21＝gk（1＋α1（kh）2＋β1（kh）4）；A22＝－ω（1＋（α1＋1／3）（kh）2＋（β1 ＋α1／3－1／45）（kh）4）；A23＝－nrω（kh）2；A31＝gk（1＋nβr（kh）2／cr）；A32＝－ω（kh）2／2；A33＝－ω（cr＋iα／ω）（1＋r2（kh）2／3）＋nr（kh）2（1＋β＋iβα／（ωcr））.

只有滿足det（A）＝0，方程（9）才有非零解，進(jìn)一步將其與Gu等的色散關(guān)系表達(dá)式對(duì)比［10］：

在保持n＝0.5不變的情況下，圖1（a）、（b）、（c）分別給出當(dāng)s＝nω／a＝0.1，r＝0.5，2.0，5.0時(shí)本文方程解（利用det（A）＝0給出的計(jì)算結(jié)果）與理論解析解（式（11）計(jì)算的結(jié)果）的對(duì)比.由圖1可見(jiàn)：

（1）在不同滲透厚度比情況下，方程的相速度與解析解吻合程度很高，這說(shuō)明方程具有較好的色散性；對(duì)比圖1（a）、（b）、（c）發(fā)現(xiàn)，當(dāng)量綱一水深h／L0＝1.0 時(shí)，隨著r的增大，誤差也逐漸增大，最大誤差為2.5%（圖1（c））.而不存在滲透介質(zhì)時(shí)，由方程（1）和方程（2）組成的Boussinesq方程在h／L0＝1.0時(shí)的最大誤差為2%，這反映出滲透介質(zhì)厚度對(duì)相速度的影響相對(duì)較小.

圖1 模型量綱一相速度、衰減率與線性波浪理論解析解的比較Fig.1 Comparisons of dimensionless phase velocity and damping rate from the model and the linear wave theory

（2）在衰減性方面，r＝2.0時(shí)，在h／L0＜0.6范圍內(nèi)，方程的衰減率與解析解相差不大，超過(guò)這一范圍后誤差迅速增大；而在其他情況下，尤其是r＝5.0時(shí)，當(dāng)h／L0超過(guò)0.4以后，本文方程的結(jié)果與解析解存在較大誤差.因此，理論分析表明，方程的衰減率與解析解的誤差相對(duì)較大，縱觀不同滲透厚度比下的情況，其可接受的最大h／L0為0.4.本文滲透方程的最高階僅為2階，這應(yīng)該是導(dǎo)致衰減率對(duì)比效果不佳的主要原因.

3 一維數(shù)值模型

采用類(lèi)似劉忠波等的方法［1］，對(duì)一維情況下的控制方程（3）～（5）進(jìn)行數(shù)值求解.即在非交錯(cuò)網(wǎng)格上離散方程，采用高精度差分格式對(duì)方程中的時(shí)間和空間導(dǎo)數(shù)進(jìn)行近似，采用三階Adams－Bashforth預(yù)報(bào)格式、四階Adams－Moulton 校正格式進(jìn)行時(shí)間積分.這里也將方程（3）～（5）整理成文獻(xiàn)［1］中的格式，具體表達(dá)如下：

式中

式中：f（x，t）是內(nèi)部造波源項(xiàng)，具體表達(dá)形式可參見(jiàn)文獻(xiàn)［11］.

進(jìn)行時(shí)間積分時(shí)，預(yù)報(bào)格式為

校正格式為

求解動(dòng)量方程時(shí)，采用五對(duì)角寬帶解法［1］.而關(guān)于時(shí)間t求導(dǎo)的項(xiàng)、為增加程序穩(wěn)定性而采用的濾波技術(shù)以及邊界消波技術(shù)，均采用文獻(xiàn)［2］給出的方法.當(dāng)兩次計(jì)算的值（3個(gè)變量）相對(duì)誤差均小于0.000 1時(shí)，當(dāng)前迭代結(jié)束，否則，重新利用式（18）～（20）進(jìn)行校正計(jì)算，主要計(jì)算流程見(jiàn)圖2.

圖2 計(jì)算流程示意圖Fig.2 The schematic diagram of calculation flow chart

4 數(shù)值驗(yàn)證

Hsiao等進(jìn)行了波浪在滲透潛堤上傳播變形的實(shí)驗(yàn)［5］.實(shí)驗(yàn)中，前坡坡度為1∶8.89，而后坡坡度為1∶5.93，最大水深為0.175m，三角形潛堤上最小水深為0.04m，潛堤由直徑約為1.9cm的石子組成，石子間孔隙率n＝0.42；實(shí)驗(yàn)中的波周期為1s，波高為2.7cm.計(jì)算中，采用Hsiao給出的參數(shù)值：α1＝5.83s－1，α2＝58.85m，時(shí)間步長(zhǎng)0.01s，空間步長(zhǎng)為0.02m.

數(shù)值結(jié)果與實(shí)驗(yàn)結(jié)果的對(duì)比見(jiàn)圖3，由圖可見(jiàn)，當(dāng)不考慮滲透作用時(shí)（文獻(xiàn)［1］數(shù)值模擬結(jié)果），在x＝10.8 m 處及以后，出現(xiàn)較大的誤差；隨著滲透作用的進(jìn)一步積累，數(shù)值結(jié)果與實(shí)驗(yàn)結(jié)果的相位差越來(lái)越大，而且波幅也存在較大差異.采用本文模型模擬時(shí)，在大多位置處，二者在相位和幅值上均較為吻合.而在個(gè)別位置處，如x＝11.4m 處相位開(kāi)始出現(xiàn)一定的誤差，且幅值上也存在一定差異，所計(jì)算的結(jié)果比實(shí)驗(yàn)結(jié)果偏小.總之，對(duì)比計(jì)算結(jié)果可知，本文模型計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的吻合程度更佳，這說(shuō)明多引入一個(gè)動(dòng)量方程來(lái)適應(yīng)滲透介質(zhì)的影響的做法是有效的.

圖3 數(shù)值計(jì)算的波面升高與實(shí)驗(yàn)結(jié)果的比較Fig.3 Comparisons of the computed surface elevation and experimental results

為進(jìn)一步驗(yàn)證模型的適用性，以Sawaragi等的實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證本文模型［12］.實(shí)驗(yàn)中的水槽長(zhǎng)30 m，寬0.7 m，高0.9 m.水槽中鋪設(shè)一段長(zhǎng)3.5 m，厚0.15m 的平底粗顆粒海床，泥沙顆粒粒徑采用了1.8cm 和3.07cm 兩種，海床上的水深為0.15m.本文選取入射波高Hi為3.58cm，周期T＝1.0s，孔 隙 率n為0.4 進(jìn) 行 數(shù) 值 模 擬.Karunarathna等曾采用該實(shí)驗(yàn)數(shù)據(jù)對(duì)Navier－Stokes數(shù)學(xué)模型進(jìn)行驗(yàn)證［13］，實(shí)驗(yàn)數(shù)據(jù)取自Karunarathna和Lin的文獻(xiàn).數(shù)值模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的比較見(jiàn)圖4，由圖可見(jiàn)，數(shù)值結(jié)果（波高）沿波浪傳播方向存在一定的震蕩，但整體來(lái)看，與實(shí)驗(yàn)數(shù)據(jù)的吻合程度良好，這反映出本文模型也能較好地模擬波浪在滲透海床上的傳播變形.

圖4 滲透海床上波浪衰減的數(shù)值結(jié)果與實(shí)驗(yàn)結(jié)果的比較（T＝1.0s）Fig.4 Comparison of numerical results against experimental data for the wave damping over porous seabed（T ＝1.0s）

5 結(jié) 論

在一組具有四階色散精度的Boussinesq 方程基礎(chǔ)上，直接耦合Cruz等導(dǎo)出的滲透介質(zhì)中流體的動(dòng)量方程，構(gòu)建了一組新的方程.理論分析了方程的相速度和衰減率，并與解析解進(jìn)行了對(duì)比，結(jié)果表明該方程可期望用于模擬波浪在滲透地形上的傳播變形；進(jìn)一步在非交錯(cuò)網(wǎng)格下對(duì)方程的一維形式進(jìn)行離散，建立了基于有限差分法的數(shù)值模型，利用數(shù)值模型模擬了波浪在滲透地形上的傳播變形，并將計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比，二者的吻合程度較高.這也從側(cè)面反映出本文的改進(jìn)方法是有效的，為非滲透地形下的其他Boussinesq類(lèi)方程擴(kuò)展為適用于滲透海床情況提供了重要參考.

［1］ 劉忠波，鄒志利，孫昭晨.加強(qiáng)的適合復(fù)雜地形的水波方程及其一維數(shù)值模型驗(yàn)證［J］.海洋學(xué)報(bào)，2008，30（3）：117－125.LIU Zhong－bo，ZOU Zhi－li，SUN Zhao－chen.Enhanced Boussinesq equations for rapidly varying topographies and their one－dimensional numerical validation ［J］.Acta Oceanologica Sinica，2008，30（3）：117－125.（in Chinese）

［2］ Kirby J T.Boussinesq models and applications to nearshore wave propagation，surf zone processes and wave－induced currents ［M］／／ Advances in Coastal Engineering.New York：Elsevier Science，2002：1－41.

［3］ Madsen P A，F(xiàn)uhrman D R.High－order Boussinesq－type modeling of nonlinear wave phenomena in deep and shallow water［M］／／Advances in Numerical Simulation of Nonlinear Water Waves.Singapore：World Scientific Publishing Co.Pte.Ltd.，2010：245－285.

［4］ Cruz E C，Isobe M，Watanabe A.Boussinesq equations for wave transformation on porous beds［J］.Coastal Engineering，1997，30：125－156.

［5］ Hsiao S，Liu P L F，Chen Y，et al.Nonlinear water waves propagating over a permeable bed［J］.Proceedings of Royal Society，London A，2002，458：1291－1322.

［6］ Chen Q.Fully nonlinear Boussinesq－type equations for waves and currents over porous beds ［J］.Journal of Engineering Mechanics，2006，132（2）：220－230.

［7］ Peregrine D H.Long waves on a beach［J］.Journal of Fluid Mechanics，1967，27（4）：815－827.

［8］ 劉忠波，孫昭晨.波浪在滲透海床上傳播的數(shù)學(xué)模型［J］.中國(guó)科技論文在線，2011，6（5）：374－379.LIU Zhong－bo，SUN Zhao－chen.Wave propagating model over a porous seabed ［J］.China Science Paper，2011，6（5）：374－379.（in Chinese）

［9］ Avgeris I，Karambas T V，Prinos P.Boussinesq modeling of wave interaction with porous submerged breakwaters ［C］ ／／ Proceedings of 29th International Conference on Coastal Engineering.Lisbon：ASCE，2004：604－616.

［10］ Gu Z，Wang H.Gravity waves over porous bottoms［J］.Coastal Engineering，1991，15：497－524.

［11］ Gobbi M F，Kirby J T.Wave evolution over submerged sills：Tests of a high－order Boussinesq model［J］.Coastal Engineering，1999，37：57－96.

［12］ Sawaragi T，Deguchi I.Waves on permeable layers［C］／／Proceedings of 23rd Conference on Coastal Engineering.Venice：ASCE，1992：1531－1544.

［13］ Karunarathna S，Lin P.Numerical simulation of wave damping over porous seabeds［J］.Coastal Engineering，2006，53：845－855.