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

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

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

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

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

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

1 數(shù)學(xué)模型

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

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

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

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

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

1.2 方程空間離散

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

1.3 數(shù)值通量計(jì)算

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

第一步:通量計(jì)算

若k=K，流程結(jié)束。

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

第三步:返回第一步。

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

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

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

1.4 波浪破碎和水-陸動(dòng)邊界處理

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

1.5 時(shí)間積分以及邊界條件

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

2 模型驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 潛礁地形上孤立波傳播

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

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

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

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

3 結(jié)語

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

［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.