房克照，劉忠波，唐軍，鄒志利，尹繼偉

(1.大連理工大學(xué)海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧大連116024;2.國(guó)家海洋局海洋環(huán)境檢測(cè)中心，遼寧大連116023;3.大連海事大學(xué)交通運(yùn)輸管理學(xué)院，遼寧大連116026;4.長(zhǎng)沙理工大學(xué)水沙科學(xué)與水災(zāi)害防治湖南省重點(diǎn)實(shí)驗(yàn)室，湖南長(zhǎng)沙410076)

數(shù)值模擬研究海嘯波的傳播以及爬高等問題具有十分重要的意義。當(dāng)近岸潛礁存在時(shí)，由于其前坡較為陡峭且潛礁內(nèi)平臺(tái)水深較小，海嘯波的傳播形態(tài)和能量組成將產(chǎn)生顯著變化，較一般開闊海岸上更為復(fù)雜［1］。數(shù)值模擬時(shí)，不僅要求模型能夠描述非線性海嘯波的傳播過程，還要具備處理波浪破碎和海岸動(dòng)邊界問題的能力。此外，海嘯波在近岸潛礁上傳播容易形成水躍，具有極其陡峭的波面形態(tài)，要求模型具有較強(qiáng)的間斷捕捉能力［1-2］。

通常采用非線性淺水方程結(jié)合孤立波數(shù)值研究海嘯波的傳播、破碎和爬高過程［3-4］。非線性淺水方程是典型的雙曲型模型，已有大量的高精度數(shù)值求解格式［5］。不但擁有高分辨率特征，同時(shí)能較方便地處理波浪破碎和海岸動(dòng)邊界問題，基于非線性淺水方程和利用高精度數(shù)值方法建立的海嘯波傳播數(shù)學(xué)模型取得了較大的成功［3-4］。但該方程存在重大缺陷:缺乏色散性，使得其僅能應(yīng)用于模擬近岸海嘯波傳播，大大限制了其應(yīng)用范圍。而近幾十年來(lái)發(fā)展得的Boussinesq類水波方程［6］，雖然精度各異，但其既包含色散性又包含非線性，具備模擬海嘯波傳播的潛質(zhì)。但是這類模型對(duì)于波浪破碎和海岸動(dòng)邊界的處理都是近似的，同時(shí)引入大量可調(diào)參數(shù)，增加了數(shù)值計(jì)算結(jié)果的不確定性［6］。此外，這類模型多應(yīng)用有限差分方法求解，計(jì)算過程中容易產(chǎn)生高頻數(shù)值震蕩。為克服這一問題，通常采用數(shù)值光滑技術(shù)［6-9］。即便較低頻率地使用光滑函數(shù)，也難以估量其對(duì)真正物理現(xiàn)象的影響程度。

結(jié)合上述2類模型的優(yōu)缺點(diǎn)，發(fā)展一種適用于海嘯波模擬且同時(shí)具備2種模型優(yōu)勢(shì)的數(shù)值模式非常有吸引力。近幾年，已經(jīng)有一些學(xué)者開展了這方面研究［10-14］。這些模型均采用早期發(fā)展的Boussinesq類方程［15-16］作為控制方程，采用有限體積和有限差分方法進(jìn)行數(shù)值求解。其中，有限體積方法在計(jì)算界面通量時(shí)均采用一階精度的通量函數(shù)。本文在Kim等拓展的一組適合陡坡地形上的Boussinesq方程［17］基礎(chǔ)上，采用具有TVD性質(zhì)的二階通量函數(shù)算法，建立基于有限體積和有限差分方法的數(shù)值模型，模擬孤立波在潛礁上的傳播。

1 控制方程

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

文獻(xiàn)［17］給出了一組能考慮快速變化水底地形的Boussinesq水波方程，其一維形式［17］如下:

其中，

式中:h為靜水深，η為波面升高，d=h+η為當(dāng)?shù)厮?，q=du為單位流量(u為水深積分平均流速);下標(biāo)x和t分別表示變量對(duì)空間和時(shí)間的導(dǎo)數(shù);B為色散性參數(shù)，取值1/15時(shí)，控制方程色散性為精確色散關(guān)系的Padé［2，2］階近似，適用于深水波浪傳播問題。在上述方程中忽略式(3)中后2項(xiàng)中與底坡坡度有關(guān)的項(xiàng)，則可退化為文獻(xiàn)［15］給出的方程，高階地形導(dǎo)數(shù)項(xiàng)的存在可以有效提高模型在快速變化地形條件下的精度［17］。

1.2 守恒形式的控制方程

為便于采用有限體積法求解控制方程，將其寫為如下守恒形式:

為取得靜水平衡解，已將式(2)中dghx改寫為［g(h2+2hη)/2］x-ghηx。Sd為色散項(xiàng)，R=-fu|u|代表水底摩擦，f為底摩擦系數(shù)，取值在后文給出。

2 數(shù)值求解方法

2.1 方程離散

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

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

圖1給出了WAF方法的示意圖，其是一種具有二階精度的Godunov類方法［5］，核心思想在于采用nΔt/2時(shí)刻計(jì)算界面數(shù)值通量，具體表達(dá)形式如下

其對(duì)應(yīng)的離散形式為

式中，SR和SL分別為左右特征波速度，按照下式計(jì)算［5］

式中，“*”號(hào)表示中間狀態(tài)變量，按如下計(jì)算［5］

式(10)中wk為各通量分量的權(quán)系數(shù)［5］，定義為

式中，ck第k個(gè)特征波對(duì)應(yīng)的CFL數(shù)，由下式給出［5］

式中，Sk為第k個(gè)特征波的波速，見圖1、2。

圖1 WAF方法示意圖Fig.1 The sketch of WAF method

圖2 局部黎曼問題示意圖Fig.2 The sketch of local Riemann problem

式中，為限制因子，本文取Minmod限制因子［5］。

為提高格式精度，對(duì)局部黎曼問題，左右狀態(tài)變量通過四階狀態(tài)插值方法(monotone upstream-centered schemes for conservation laws，MUSCL)［8］進(jìn)行。

2.3 離散方程的時(shí)間積分

時(shí)間積分由 MUSCL-Hancock方式［14］進(jìn)行，該方法為預(yù)估校正兩步法。其中，預(yù)測(cè)步為

校正步為利用預(yù)測(cè)部給出的數(shù)值解進(jìn)行WAF數(shù)值通量計(jì)算

時(shí)間步長(zhǎng)的選取則滿足如下CFL條件

式中，ν取 0.35。

根據(jù)人才培養(yǎng)方案及課程標(biāo)準(zhǔn)，希望學(xué)生通過本課程學(xué)習(xí)，能夠正確認(rèn)識(shí)艦艇環(huán)境因素與人體健康的關(guān)系，進(jìn)一步強(qiáng)化預(yù)防為主的理念，了解海軍衛(wèi)生學(xué)在海軍建設(shè)中的重要作用，樹立衛(wèi)生學(xué)是維護(hù)和提高部隊(duì)作戰(zhàn)（作業(yè)）能力的重要保障這一專業(yè)信念，掌握發(fā)現(xiàn)和解決部隊(duì)平戰(zhàn)時(shí)各種衛(wèi)生學(xué)問題必需的基礎(chǔ)知識(shí)及技能，為今后從事軍事醫(yī)學(xué)各個(gè)領(lǐng)域的工作，尤其是開展海軍衛(wèi)生學(xué)工作打下良好基礎(chǔ)。

2.4 邊界條件

本文涉及到固壁(完全反射邊界)、吸收邊界、海岸動(dòng)邊界。

對(duì)于固體邊界，采用外設(shè)3層虛擬網(wǎng)格的方法處理(也考慮到四階MUSCL方法需要左右3個(gè)點(diǎn)網(wǎng)格點(diǎn)值)，即

吸收邊界則通過設(shè)置海綿層實(shí)現(xiàn)，即計(jì)算結(jié)果乘一光滑函數(shù):

式中，xs為海綿層厚度。

當(dāng)涉及到海岸動(dòng)邊界問題時(shí)，對(duì)式(12)、(13)中特征速度進(jìn)行修正［2，8］:

其中，干濕網(wǎng)格的分界定義為水深0.002 m。

根據(jù)文獻(xiàn)［11］建議，當(dāng)波面升高同水深比達(dá)到0.8時(shí)，認(rèn)為波浪發(fā)生破碎，控制方程(1)～(3)中的Boussinesq非線性項(xiàng)和色散性不參與計(jì)算，方程退化為淺水非線性方程，波浪破碎處理為激波。

3 計(jì)算結(jié)果討論和分析

本節(jié)針對(duì)所建立數(shù)學(xué)模型進(jìn)行驗(yàn)證和應(yīng)用。數(shù)值算例包括常水深水槽中強(qiáng)非線性孤立波的傳播和3種代表性潛礁地形上孤立波的傳播。數(shù)值結(jié)果將與解析解或?qū)嶒?yàn)數(shù)據(jù)進(jìn)行分別對(duì)比。

3.1 孤立波在常水深水槽中的傳播

數(shù)值模擬孤立波在常水深水槽中的傳播是檢驗(yàn)Boussinesq類數(shù)學(xué)模型的經(jīng)典算例［7］。孤立波在傳播過程中保持形狀和速度不變，這是非線性和色散性相互制約并達(dá)到平衡的結(jié)果。雖然Boussinesq類模型自身既具有色散性又具有非線性，但若數(shù)值格式不合適則會(huì)引入額外偽色散或耗散等，從而難以準(zhǔn)確模擬孤立波傳播過程。

圖3 不同時(shí)刻孤立波形狀數(shù)值計(jì)算結(jié)果Fig.3 The computed solitary wave profile at different moments

圖4 t=60 s時(shí)數(shù)值計(jì)算結(jié)果同解析解的對(duì)比Fig.4 Comparisons of surface elevation between computed results and analytic solution at t=60 s

數(shù)值模擬時(shí)，水槽長(zhǎng)度375 m，水深h=1.0 m，孤立波波高A=0.6 m，A/h=0.6，屬于強(qiáng)非線性波浪。空間網(wǎng)格尺寸Δx=0.05 m，計(jì)算時(shí)在計(jì)算域內(nèi)給定孤立波解析解(位于x=50 m處)。圖3給出了8個(gè)不同時(shí)刻，孤立波波面的計(jì)算結(jié)果，可以看出孤立波波形在長(zhǎng)時(shí)間內(nèi)保持不變。圖4中給出t=70 s時(shí)，數(shù)值解同解析解的對(duì)比，兩者幾乎重合，這表明孤立波傳播速度模擬準(zhǔn)確。本算例表明，所建立數(shù)值格式有效，可用于模擬孤立波傳播問題.

3.2 孤立波在不同潛礁地形上的傳播

為了研究海嘯波在潛礁地形上的傳播特性，Roeber進(jìn)行了孤立波在不同潛礁地形上的物理模型實(shí)驗(yàn)［2］，本節(jié)將針對(duì)其中的3個(gè)典型工況進(jìn)行模擬，各工況設(shè)置和網(wǎng)格尺寸在表1中給出。計(jì)算時(shí)，底摩擦系數(shù)取0.007 5，左側(cè)設(shè)置海綿層5 m，右側(cè)為完全反射邊界。

表1 各工況設(shè)置以及網(wǎng)格尺寸Table 1 Case setting and grid size in simulation

為定量考察計(jì)算結(jié)果和實(shí)驗(yàn)數(shù)據(jù)的吻合程度，引入 Wilmott因子［18］:

式中，下標(biāo)v為要考察的變量，Y(j)和y(j)分別為計(jì)算和實(shí)測(cè)數(shù)據(jù)，為實(shí)測(cè)數(shù)據(jù)平均值。當(dāng)Iv=1時(shí)表示兩者完全吻合，Iv=0時(shí)表示兩者不吻合。

工況1的計(jì)算結(jié)果如圖5。該工況中，潛礁初始處于無(wú)水狀態(tài)，因此涉及到動(dòng)邊界的處理，對(duì)模型處理干濕邊界的能力是一個(gè)考驗(yàn)。t＇=52.79，波浪開始在初始干潛礁上傳播，在x=22 m處，流體經(jīng)歷從亞臨界流速轉(zhuǎn)化為超臨界流速狀態(tài)［2］。t＇=54.79 開始，波浪開始崩塌，波前迅速超前傳播，同時(shí)一部分反射波浪往相反方向傳播，t＇=65.44時(shí)刻，由于兩部分水體的相反運(yùn)動(dòng)，部分潛礁前坡露出水面?？傮w而言，數(shù)值結(jié)果同實(shí)驗(yàn)數(shù)據(jù)吻合較好，驗(yàn)證了模型的正確性，尤其是反映出本文模型具有處理干濕動(dòng)邊界問題的能力。同時(shí)經(jīng)計(jì)算，其Wilmott因子0.975，這表明數(shù)值結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合度較高。

圖5 工況1計(jì)算波面同實(shí)驗(yàn)數(shù)據(jù)的對(duì)比Fig.5 Comparisons of computed and measured surface elevations for case 1

圖6給出了工況2的計(jì)算結(jié)果。與工況1不同，該算例潛礁一直處于淹沒狀態(tài)。t＇=67.1開始，出現(xiàn)激波，即波浪破碎現(xiàn)象。激波形成后一直保持尖銳的形狀向前傳播，至固壁邊界后被完全反射后仍然保持尖銳形狀向相反方向傳播。t＇=108.5時(shí)刻，反射波傳入深水，色散性開始發(fā)揮作用，波浪分散成不同形態(tài)的短波。數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合較好，表明模型具有良好的處理激波的能力。該組模擬對(duì)應(yīng)的Wilmott因子為0.965，說明數(shù)值結(jié)果同實(shí)驗(yàn)數(shù)據(jù)吻合度較高。

工況3在初始時(shí)刻潛礁峰露出水面0.06 m，而潛礁平臺(tái)淹沒于水下0.14 m，地形設(shè)置類似自然界中帶瀉湖的潛礁，非常具有代表性。計(jì)算結(jié)果如圖7。可見，孤立波在潛礁前坡傳播過程中，波前變陡，破碎發(fā)生形成激波。t＇=59.31時(shí)刻起，波浪開始跨越潛礁峰傳播至瀉湖區(qū)域。在潛礁峰后激起水躍，并保持尖銳形狀向前傳播。被完全反射后反向傳播，至深水區(qū)域由于色散性作用加強(qiáng)也出現(xiàn)諸多短波。數(shù)值結(jié)果同實(shí)驗(yàn)數(shù)據(jù)吻合良好，表明模型也適用于模擬這種特殊地形條件下孤立波傳播和破碎現(xiàn)象。計(jì)算得到Wilmott因子為0.942，再次表明模型具有較高的精度。

圖6 工況2計(jì)算波面同實(shí)驗(yàn)數(shù)據(jù)的對(duì)比Fig.6 Comparisons of computed and measured surface elevations for case 2

圖7 工況3計(jì)算波面同實(shí)驗(yàn)數(shù)據(jù)的對(duì)比Fig.7 Comparisons of computed and measured surface elevations for case 3

在圖5～7中也給出了 Roeber的計(jì)算結(jié)果［2］，本文計(jì)算結(jié)果和其較為接近，表明本文模型的數(shù)值離散和求解過程正確，也表明色散性Boussinesq類水波方程適用于描述潛礁地形上破碎波傳播問題。二者的差別主要由以下2個(gè)因素造成:1)模型采用的控制方程不同，Roeber模型建立在Nowgu方程［19］基礎(chǔ)上，本文采用Kim給出的方程［17］;2)數(shù)值格式不同，Roeber模型采用HLL方法［5］計(jì)算通量，采用ABM方法進(jìn)行時(shí)間積分，而本文采用WAF方法計(jì)算通量，采用MUSCL-Hancock方式進(jìn)行時(shí)間積分。

4 結(jié)論

本文建立了求解Boussinesq水波方程的數(shù)值模式，用于模擬具有孤立波形狀的海嘯波在潛礁地形上的傳播、爬坡和破碎過程.得到如下結(jié)論:

1)較傳統(tǒng)有限差分格式而言，采用有限體積和有限差分混合方法建立的模型，具備數(shù)值穩(wěn)定、捕捉激波、動(dòng)邊界處理簡(jiǎn)單以及可調(diào)參數(shù)少的優(yōu)點(diǎn)，更適用于模擬孤立波在潛礁地形上傳播的復(fù)雜過程。

2)數(shù)值解與解析解和實(shí)驗(yàn)數(shù)據(jù)吻合較好，表明所建立模型具有較高的精度，同時(shí)非常適用于潛礁地形上海嘯波的數(shù)值模擬研究。

3)所建立數(shù)值方法也為求解其他Boussinesq方程的數(shù)值提供了參考。

［1］YAO Y，HUANG Z H，MONISMITH S G，et al.DH Boussinesq modeling of wave transformation over fringing reefs［J］.Ocean Engineering，2012，47:30-42.

［2］ROBER V.Boussinesq-type model for nearshore wave processes［D］.Honolulu:University of Hawaii at Manoa，2010:47-122.

［3］HUBBARD M E，DODD N.A 2D numerical model of wave runup and overtopping［J］.Coastal Engineering，2002，47:1-26.

［4］LIANG Q，DU G Z.Flood inundation modeling with an adaptive quadtree grid shallow water equation solver［J］.Journal of Hydraulic Engineering，2008，134(1):1603-1610.

［5］TORO E F.Riemann solvers and numerical methods for fluid dynamics:a practical introduction［M］.New York:Springer Press，2009:413-578.

［6］KIRBY J T.Boussinesq models and applications to nearshore wave propagation，surfzone processes and wave-induced currents［M］.New York:Elsevier Science，2002:1-41.

［7］WEI G，KIRBY J T，GRILLI S T，et al.Full nonlinear Boussinesq model for surface waves.I:Highly nonlinear，unsteady waves［J］.Journal of Fluid Mechanics，1995(294):71-92.

［8］SHI F Y，KIRBY J T，HARRIS J C，et al.A higher-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation［J］.Ocean Modelling，2012，43/44:36-51.

［9］LYNETT J.Nearshore wave modeling with higher-order Boussinesq-type equations［J］.Journal of Waterway，Port，Coastal，and Ocean Engineering，2006，132:348-357.

［10］ERDURAN K S.Further application of hybrid solution to another form of Boussinesq equations and comparisons［J］.International Journal for Numerical Methods in Fluids，2007，53:827-849.

［11］TONELLI M，PETTI M.Hybrid finite volume-finite difference scheme for 2DH improved Boussinesq equations［J］.Coastal Engineering，2009，56:609-620.

［12］ORSZAGHOVA J，BORTHWICK AGL，TAYLOR P H，et al.From the paddle to the beach-A Boussinesq shallow water numerical wave tank based on Madsen and S?rensen’s equations［J］.Journal of Computational Physics，2011，231:328-344.

［13］CIENFUEGOS R.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.

［14］SHIACH J B，MINGHAM C G.A temporally second-order accurate Godunov-type scheme for solving the extended Boussinesq equations［J］.Coastal Engineering，2009，56:32-45.

［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］NWOGU O.An alternative form of the Boussinesq equations for near shore wave propagation ［J］.Journal of Waterway，Port，Coastal and Ocean Engineering，1993，119(6):618-638.

［17］KIM G，LEE C，SUH K D.Extended Boussinesq equations for rapidly varying topography［J］.Ocean Engineering，2009，36:842-851.

［18］WILMOTT C J.On the validation of models［J］.Physical Geograhpy，1981，2:184-194.

［19］NWOGU O.An alternative form of the Boussinesq equations for nearshore wave propagation［J］.Journal of Waterway，Port，Coastal，and Ocean Engineering，1993，119(6):618-638.