Boussinesq水波方程新型數值解法

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

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

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

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

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

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

1 數學模型

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

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

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

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

式中:f為底摩擦系數。

1.2 方程空間離散

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

1.3 數值通量計算

高分辨率數值通量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格式具有簡潔、直觀的表達形式，不需要求解局部黎曼解的特征波狀態。而傳統的HLL格式［1-4，6-9］需要判斷黎曼解特征波狀態，需要計算特征波波速，并以此進行通量計算，計算繁瑣。因此，MUSTA格式在計算效率和程序編制方面均要比HLL格式具有優勢。

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

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

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

1.5 時間積分以及邊界條件

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

2 模型驗證

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

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

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

針對這一問題，文獻［8］給出了解析解。t=100 s時，數值解同解析解的對比在圖3給出，可見數值解同解析解高度吻合。此外，數值計算的傳播速度為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格式的區別，引入以下誤差判斷標準［13］:

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

表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 規則波在潛堤地形上的傳播

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

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

2.3 潛礁地形上孤立波傳播

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

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

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

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

3 結語

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