Boussinesq水波方程新型數(shù)值解法

房克照，孫家文，劉忠波，尹 晶，張 哲

(1.大連理工大學海岸和近海工程國家重點實驗室，遼寧大連 116024;2.國家海洋環(huán)境監(jiān)測中心，遼寧大連 116023)

由于在模型穩(wěn)定性、波浪破碎處理和海岸動邊界問題方面較傳統(tǒng)有限差分方法具有優(yōu)勢，近十年來，將有限體積方法應用于Boussinesq類水波方程的數(shù)值求解成為研究熱點，形成了一類有限體積和有限差分方法混合求解的數(shù)值格式［1-9］。

建立該類混合數(shù)值格式的主要思路如下，將Boussinesq類方程寫為守恒格式，對通量項采用有限體積方法計算，剩余項采用有限差分方法求解。當波浪破碎在水-陸交界附近，色散性和高階非線性對流體運動貢獻較小，不參與計算，則Boussinesq方程退化為完全非線性淺水方程。完全非線性淺水方程屬于典型雙曲守恒律方程，有許多基于有限體積方法的高分辨率數(shù)值格式［10］(多以黎曼間斷解問題為理論基礎構筑)，因此其足以勝任破碎波浪(視為水躍)和水-陸交界(視為間斷解)的數(shù)值計算。

由于Boussinesq方程中大多數(shù)項仍通過有限差分方法求解，有限體積方法的數(shù)值實現(xiàn)重點體現(xiàn)在網(wǎng)格界面處單調(diào)性數(shù)值通量的計算。通常，采用中心格式和迎風格式進行［10-12］。中心格式不需要詳細了解特征波傳播狀態(tài)，構造簡單，但精度低。迎風格式(如HLL格式和Roe格式)精度高，但需要精確或者近似求解黎曼問題，計算量大，構造過程相對繁瑣?，F(xiàn)有混合格式大都采用HLL格式進行數(shù)值通量計算［1-4，6-9］。也有部分學者采用Roe格式計算［5］，但其每個時間步都涉及到特征矩陣的求解和分裂計算，十分繁瑣。HLL和Roe格式均建立于詳細的黎曼解狀態(tài)基礎上，并且均存在一些不足，譬如HLL格式在水平二維問題時需要擴展至HLLC格式以彌補不能考慮接觸波的不足，Roe格式某些情況下不滿足熵增條件，對于實際應用而言，較準確地求解黎曼解不是一件易事［10］。Toro等近年來發(fā)展了MUSTA格式［11-12］，與傳統(tǒng)中心格式和迎風格式不同，其通過控制方程直接構造數(shù)值通量，是一種既具有中心格式的簡單性同時具備迎風格式精確性的新型數(shù)值格式，其精確性和有效性在求解非線性淺水方程時得到了驗證［11-13］。

將MUSTA格式應用于Boussinesq類水波方程的求解，建立有限體積/有限差分混合求解數(shù)值格式。除對所建立模型進行驗證外，重點比較MUSTA和HLL兩種格式在計算精度、計算效率以及程序編制等方面的表現(xiàn)。

1 數(shù)學模型

1.1 適用于快速變化地形的Boussinesq水波方程

以波面η和通量q表達的二階Boussinesq方程［14］為:

式中:h為靜水深;d=h+η為當?shù)厮?g為重力加速度;q=du為斷面流量(u為水深平均速度)，下標x和t分別表示變量對空間和時間的偏導數(shù);B=1/15為自由參數(shù)，通過優(yōu)化方程的色散性得到［15］。上述方程為Madsen和Sorensen方程［15］的擴展，更適用于在地形快速變化水域中應用。忽略式(2)中色散項ψx，方程退化為完全非線性淺水方程。

考慮海底摩擦，并將上述方程寫成守恒形式如下:

式中:f為底摩擦系數(shù)。

1.2 方程空間離散

將計算域在空間、時間上做如下離散:xi=i(i=1，…，N)，tn=nΔt，其中Δx、Δt分別為空間、時間網(wǎng)格步長。在有限體積［xi-1/2，xi+1/2］×［tn，tn+1］內(nèi)對控制方程(4)進行積分并應用格林定理，可得

1.3 數(shù)值通量計算

高分辨率數(shù)值通量Fi+1/2的計算均通過求解黎曼解問題進行，即求解初始值問題:

第一步:通量計算

若k=K，流程結束。

第二步:通過控制方程求近似黎曼解

第三步:返回第一步。

圖1 MUSTA格式構造示意圖(上標表示第k步值)Fig.1 The sketch of MUSTA scheme(subscript k denotes the value after kth stage)

文中，k=1，2，3時對應的格式分別簡稱為MUSTA-1，MUSTA-2和MUSTA-3。MUSTA格式具有簡潔、直觀的表達形式，不需要求解局部黎曼解的特征波狀態(tài)。而傳統(tǒng)的HLL格式［1-4，6-9］需要判斷黎曼解特征波狀態(tài)，需要計算特征波波速，并以此進行通量計算，計算繁瑣。因此，MUSTA格式在計算效率和程序編制方面均要比HLL格式具有優(yōu)勢。

上述MUSTA格式僅為一階精度，為提高精度，采用四階精度的狀態(tài)插值方法［17］對界面左右變量進行重構，利用重構后的變量進行上述計算。

1.4 波浪破碎和水-陸動邊界處理

在混合格式中，波浪破碎發(fā)生時，Boussinesq方程退化為完全非線性淺水方程，處理為水躍。根據(jù)建議［1］，當波面升高同水深比達到0.8時，認為波浪發(fā)生破碎，控制方程(1)～(3)中的Boussinesq高階非線性項和色散性不參與計算，方程退化為淺水非線性方程。對于水-陸動邊界，采用簡單有效的薄層水體法處理［13］，即定義一極小水深，本文取0.001 m，若網(wǎng)格水深大于該值視為水域;若網(wǎng)格水深小于該值，視為干網(wǎng)格，水深強制賦值0.001 m，流速為零。

1.5 時間積分以及邊界條件

時間積分通過具有TVD性質(zhì)的三階龍格-庫塔方法進行［1］。計算域兩端設置為固壁邊界，模型采用在質(zhì)量方程中增加源項的方法產(chǎn)生波浪，同時視需要在計算域末端設置海綿吸收層吸收波浪［1］。

2 模型驗證

本節(jié)將通過幾個典型算例的數(shù)值模擬，對所建立模型進行驗證。除采用MUSTA格式進行計算外，在保持其他計算條件一致情況下，也采用HLL格式進行了計算，以便兩種格式的比較。

2.1 孤立波在等水深水槽中的傳播

孤立波是波浪色散性和非線性平衡制約的典型代表，其在常水深水槽中長距離傳播的數(shù)值模擬是檢驗混合類格式的有力工具。計算時，計算域長度500 m，水深1.0 m。在計算域內(nèi)給出孤立波作為初始條件，波幅a=0.6 m，幅值中心位于x=60 m處，網(wǎng)格尺寸Δx=0.10 m，計算域兩端設置5 m長海綿層。采用MUSTA-1格式的計算結果在圖2中給出。可見，經(jīng)過長距離的傳播，波形和幅值保持不變，表明所建立的數(shù)值格式?jīng)]有引入數(shù)值耗散以及偽色散。

針對這一問題，文獻［8］給出了解析解。t=100 s時，數(shù)值解同解析解的對比在圖3給出，可見數(shù)值解同解析解高度吻合。此外，數(shù)值計算的傳播速度為4.038 m/s，同解析解4.037 3 m/s吻合。

圖2 孤立波傳播過程中不同時刻波面升高Fig.2 The surface elevations at different moments during the solitary wave propagation

圖3 t=100 s時計算波面同解析解的比較Fig.3 Comparison of computed and analytical surface profiles at t=100 s

為定量考證MUSTA格式同HLL格式的區(qū)別，引入以下誤差判斷標準［13］:

式中:Y代表要比較的變量，下標i為網(wǎng)格標記，上標sim和exact分別為數(shù)值解和解析解。計算結果在表1中給出，同時給出各種格式的CPU耗時。對比可見，對MUSTA格式而言，隨計算步數(shù)增加，模型精度提高，但計算時間也顯著增加。

表1 不同格式誤差匯總Tab.1 A summary of error norms for various MUSTA and HLL schemes

當K=3時，計算精度的提高幾乎可以忽略。MUSTA格式計算效率均高于HLL格式。其中，同HLL格式相比，MUSTA-1精度略有降低，但計算耗時顯著減少。綜合考慮計算精度、計算效率以及程序編制難易程度，本文以下計算均采用MUSTA-1格式進行。這同其他學者［12-13］推薦在求解非線性淺水方程時使用MUSTA-1格式一致。

2.2 規(guī)則波在潛堤地形上的傳播

規(guī)則波跨越潛堤傳播是非常復雜的波浪演化過程，涉及波浪的變淺、高次諧波產(chǎn)生釋放等，是檢驗模型色散性和非線性綜合性能的有力工具。這里針對Beji和Battjes實驗［17］中工況(a)進行模擬，入射波波高H=0.02 m，周期T=2.02 s，實驗布置見圖4。計算時，空間步長Δx=0.02 m，計算域兩端設置1.5倍波長海綿吸收層。圖5給出6個浪高儀位置上計算結果同實驗數(shù)據(jù)的對比。MUSTA-1和HLL格式給出的數(shù)值結果幾乎相同。在x=13.5 m以前，數(shù)值結果同實驗數(shù)據(jù)吻合。在此之后，由于控制方程色散性的不足，模擬波面同實驗數(shù)據(jù)差別較大。經(jīng)統(tǒng)計發(fā)現(xiàn)，MUSTA-1較HLL格式CPU耗時減少約5%。

圖4 規(guī)則波跨越潛堤傳播實驗布置示意(Beji和Battjes［17］)Fig.4 The sketch of experimental setup of regular waves propagation over submerged bar(Beji and Battjes［17］)

2.3 潛礁地形上孤立波傳播

為驗證模型，針對孤立波在潛礁上傳播過程進行數(shù)值模擬，其中涉及到波浪破碎、水-陸動邊界、水躍形成傳播等復雜過程，是檢驗模型較為苛刻的算例。這里，采用Roeber實驗［2］中第三組進行模擬，該組實驗地形設置類似自然界中帶潟湖潛礁，非常具有代表性。潛礁峰露出水面0.06 m，而潛礁平臺淹沒于水下0.14 m［2］。

圖5 波面升高數(shù)值結果和實驗數(shù)據(jù)對比Fig.5 The comparison of surface elevations between numerical model and experimental data

圖6 波面升高數(shù)值結果和實驗數(shù)據(jù)對比Fig.6 The comparison of surface elevations between numerical model and experimental data

計算時，底摩擦系數(shù)取0.007 5，左側設置海綿層5 m，右側為完全反射邊界，Δx=0.10 m。入射孤立波波高0.75 m，造波板處水深2.5 m。數(shù)值模擬結果在圖6中給出?？梢?，孤立波在潛礁前坡傳播過程中，波前變陡，破碎發(fā)生形成水躍。t’=59時刻起(t’=t/(g/h0)1/2為無因次時間，h0為造波板處水深)，波浪開始跨越潛礁峰傳播至潟湖區(qū)域。在潛礁峰后激起水躍，并保持尖銳形狀向前傳播(至t’=80.11)。隨后(t’=91.26，100.61和117.11)，被完全反射后反向傳播，至深水區(qū)域由于色散性作用加強也出現(xiàn)諸多短波。數(shù)值結果同實驗數(shù)據(jù)的良好吻合，表明模型也適用于這種特殊地形條件下孤立波傳播和破碎的數(shù)值模擬。同本文其他算例一致，MUSTA-1和HLL格式給出的數(shù)值結果接近，但前者計算速度提高約6%。

3 結語

將MUSTA格式應用建立求解Boussinesq水波方程的混合數(shù)值格式。針對一維守恒格式的控制方程，采用有限體積方法求解數(shù)值通量，剩余項利用有限差分方法求解。時間積分、邊界處理等同現(xiàn)有的混合格式基本類似，但本文采用MUSTA格式計算數(shù)值通量。除進行模型驗證外，重點對MUSTA和HLL格式進行了對比研究。綜合考慮計算精度、計算效率、程序編制和實際應用，MUSTA-1格式較普遍應用的HLL格式具有優(yōu)勢，值得在Boussinesq類水波方程的求解方面進行推廣和應用。

［1］ Shi F，Kirby J T，Harris J C，et al.A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation ［J］.Ocean Modelling，2012，43-44:36-51.

［2］ Roeber V，Cheung K F.Boussinesq-type model for energetic breaking waves in fringing reef environments［J］.Coastal Engineering，2012，70:1-20.

［3］ Roeber V，Cheung K F，Kobayashi M H.Shock-capturing Boussinesq-type model for nearshore wave processes［J］.Coastal Engineering，2010，57(4):407-423.

［4］ Bradford S F，Sanders B F.Finite volume schemes for the Boussinesq equations［J］.Ocean Wave Measurement and Analysis，2001:953-962.

［5］ Erduran K S，Llic S，Kutijia V.Hybrid-finite volume finite-difference scheme for the solution of Boussinesq equations［J］.International Journal for Numerical Methods in Fluids，2005，49:1213-1232.

［6］ Cienfuegos R，Barthelemy E，Bonneton P.A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations.Part II:boundary conditions and validation ［J］.International Journal for Numerical Methods in Fluids，2007，53:1423-1455.

［7］ Shiach J B，Mingham C G.A temporally second-order accurate Godunov-type scheme for solving the extended Boussinesq equations［J］.Coastal Engineering，2009，59(1):32-45.

［8］ Orszaghova J，Borthwick A G L，Taylor P H.From the paddle to the beach-A Boussinesq shallow water numerical wave tank based on Madsen and Sorensen’s equations［J］.Journal of Computational Physics，2012，231(2):328-344.

［9］ Dutykh D，Katsaounis T，Mistotakis D.Finite volume scheme for dispersive wave propagation and runup［J］.Journal of Computational Physics，2011，230(8):3035-3061.

［10］ Toro E F.Riemann solvers and numerical methods for fluid dynamics［M］.Springer Verlag，2009.

［11］ Toro E F.Musta fluxes for systems of conservation laws［J］.Journal of Computational Physics，2006，216(2):403-429.

［12］ Toro E F.Musta:A multi-stage numerical flux［J］.Applied Numerical Mathematics，2006，56(10-11):1464-1479.

［13］ Guo W D，Lai J S，Lin G F.Finite-volume multi-stage schemes for shallow-water flow simulation［J］.International Journal for Numerical Methods in Fluids，2008，57(2):177-204.

［14］ Kim G，Lee C，Suh K.Extended Boussinesq equations for rapidly varying topography［J］.Ocean Engineering，2009，36:842-851.

［15］ Madsen P A，Sorensen O R.A new form of the Boussinesq equations with improved linear dispersion characteristics.Part 2.A slowly-varying bathymetry［J］.Coastal Engineering，1992，18:183-204.

［16］ Yamamoto S，Kano S，Daiguji H.An efficient CFD approach for simulating unsteady hypersonic shock-shock interference flows［J］.Computers＆ Fluids，1998，27(5-6):571-580.

［17］ Beji S，Battjes J A.Experimental investigation of wave propagation over a bar［J］.Coastal Engineering，1993，19(1-2):151-162.