潛礁上孤立波傳播的數值模擬

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

(1.大連理工大學海岸和近海工程國家重點實驗室，遼寧大連116024;2.國家海洋局海洋環境檢測中心，遼寧大連116023;3.大連海事大學交通運輸管理學院，遼寧大連116026;4.長沙理工大學水沙科學與水災害防治湖南省重點實驗室，湖南長沙410076)

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

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

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

1 控制方程

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

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

其中，

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

1.2 守恒形式的控制方程

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

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

2 數值求解方法

2.1 方程離散

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

2.2 數值通量計算

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

其對應的離散形式為

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

式中，“*”號表示中間狀態變量，按如下計算［5］

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

式中，ck第k個特征波對應的CFL數，由下式給出［5］

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

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

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

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

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

2.3 離散方程的時間積分

時間積分由 MUSCL-Hancock方式［14］進行，該方法為預估校正兩步法。其中，預測步為

校正步為利用預測部給出的數值解進行WAF數值通量計算

時間步長的選取則滿足如下CFL條件

式中，ν取 0.35。

根據人才培養方案及課程標準，希望學生通過本課程學習，能夠正確認識艦艇環境因素與人體健康的關系，進一步強化預防為主的理念，了解海軍衛生學在海軍建設中的重要作用，樹立衛生學是維護和提高部隊作戰（作業）能力的重要保障這一專業信念，掌握發現和解決部隊平戰時各種衛生學問題必需的基礎知識及技能，為今后從事軍事醫學各個領域的工作，尤其是開展海軍衛生學工作打下良好基礎。

2.4 邊界條件

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

對于固體邊界，采用外設3層虛擬網格的方法處理(也考慮到四階MUSCL方法需要左右3個點網格點值)，即

吸收邊界則通過設置海綿層實現，即計算結果乘一光滑函數:

式中，xs為海綿層厚度。

當涉及到海岸動邊界問題時，對式(12)、(13)中特征速度進行修正［2，8］:

其中，干濕網格的分界定義為水深0.002 m。

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

3 計算結果討論和分析

本節針對所建立數學模型進行驗證和應用。數值算例包括常水深水槽中強非線性孤立波的傳播和3種代表性潛礁地形上孤立波的傳播。數值結果將與解析解或實驗數據進行分別對比。

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

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

圖3 不同時刻孤立波形狀數值計算結果Fig.3 The computed solitary wave profile at different moments

圖4 t=60 s時數值計算結果同解析解的對比Fig.4 Comparisons of surface elevation between computed results and analytic solution at t=60 s

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

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

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

表1 各工況設置以及網格尺寸Table 1 Case setting and grid size in simulation

為定量考察計算結果和實驗數據的吻合程度，引入 Wilmott因子［18］:

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

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

圖5 工況1計算波面同實驗數據的對比Fig.5 Comparisons of computed and measured surface elevations for case 1

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

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

圖6 工況2計算波面同實驗數據的對比Fig.6 Comparisons of computed and measured surface elevations for case 2

圖7 工況3計算波面同實驗數據的對比Fig.7 Comparisons of computed and measured surface elevations for case 3

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

4 結論

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

1)較傳統有限差分格式而言，采用有限體積和有限差分混合方法建立的模型，具備數值穩定、捕捉激波、動邊界處理簡單以及可調參數少的優點，更適用于模擬孤立波在潛礁地形上傳播的復雜過程。

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

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

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