焦子峰，尹 晶,2，房克照，孫家文,2

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

一維全非線性Green-Naghdi水波方程數(shù)值模型

焦子峰1，尹 晶1,2，房克照1，孫家文1,2

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

建立了求解一維全非線性Green-Naghdi水波方程的中心有限體積/有限差分混合數(shù)值格式。采用結(jié)構(gòu)化網(wǎng)格對(duì)守恒形式的控制方程進(jìn)行離散和積分，界面數(shù)值通量采用有限體積法計(jì)算，剩余項(xiàng)則采用中心有限差分格式求解。其中，采用中心迎風(fēng)有限體積格式計(jì)算控制體界面數(shù)值通量，并結(jié)合界面變量的線性重構(gòu)方法，使其在空間上具有四階精度，通過(guò)引入靜壓重構(gòu)技術(shù)和波浪破碎指標(biāo)使模型具備處理海岸水-陸動(dòng)邊界及波浪破碎的能力。時(shí)間積分則采用具有總時(shí)間變差減小(Total Variation Diminishing, TVD)性質(zhì)的三階龍格-庫(kù)塔法進(jìn)行。應(yīng)用該模型對(duì)孤立波在常水深和斜坡海岸上的傳播過(guò)程及規(guī)則波跨越潛堤傳播的實(shí)驗(yàn)進(jìn)行了數(shù)值模型研究，數(shù)值計(jì)算同解析解及實(shí)驗(yàn)數(shù)據(jù)吻合良好。

Green-Naghdi水波方程；混合格式；中心迎風(fēng)格式；水-陸動(dòng)邊界；波浪破碎

Abstract: A 1D Green-Naghdi model based on hybrid central finite-volume/finite-difference method is presented in this paper. The convective flux in the conservative governing equations is estimated using finite-volume method while the remaining terms are discretized using finite-difference method. Central upwind finite-volume scheme, in conjunction with fourth-order MUSCL method for variables reconstruction, is used to compute the interface flux. By introducing the hydrostatic reconstruction technique and wave breaking index, the model has the ability to deal with the coastal moving shoreline and wave breaking. The third-order Runge-Kutta method with TVD property is adopted to perform time marching. The propagation of solitary wave on plane and sloping beach and the experiment of regular propagating over a submerged sill are computed with this model. The numerical results are in good agreement with analytical and experimental data.

Keywords: Green-Naghdi wave equations; hybrid scheme; central upwind scheme; moving shoreline; wave breaking

準(zhǔn)確模擬波浪地傳播變形對(duì)海岸工程和海洋工程研究具有重要的意義。盡管直接求解水波歐拉方程的數(shù)值方法在近幾年有較快地發(fā)展[1]，Boussinesq模型仍然是研究工程尺度范圍內(nèi)波浪傳播變形問(wèn)題最為廣泛的手段之一(如開(kāi)源模型FUNWAVE[2]和COULWAVE[3]及丹麥DHI公司開(kāi)發(fā)的MIKE21商業(yè)軟件波浪模塊等)。早期的Boussinesq方程通常僅具有弱非線性特征，只適用于小波幅波浪地模擬，經(jīng)過(guò)幾十年的研究，目前Boussinesq類水波方程的相關(guān)理論已達(dá)到了一個(gè)較高的水平，為強(qiáng)非線性波浪的數(shù)值計(jì)算奠定了基礎(chǔ)[4-5]。Green-Naghdi水波方程[6]是較為經(jīng)典的、具有全非線性特征的Boussinesq類水波方程，可用于研究非線性強(qiáng)和水深相對(duì)較深條件下波浪傳播問(wèn)題(僅具有弱色散性)，因此，近幾年基于Green-Naghdi水波方程展開(kāi)的數(shù)值計(jì)算層出不窮[7-8]。由于Boussinesq類方程通常含有高階導(dǎo)數(shù)項(xiàng)(包括色散項(xiàng)和非線性項(xiàng))，不便于采用諸如有限元和有限體積方法求解，一般采用有限差分方法進(jìn)行求解[9]。雖然有限差分方法簡(jiǎn)單明了、易于編程，但采用有限差分方法的數(shù)學(xué)模型易產(chǎn)生高頻數(shù)值震蕩，當(dāng)非線性較強(qiáng)或者海床地形變化較大時(shí)尤為突出，經(jīng)常導(dǎo)致計(jì)算失敗，即便采用光滑器穩(wěn)定計(jì)算仍然不能從根本上解決上述問(wèn)題[10]。且Boussinesq類模型在處理海岸水-陸動(dòng)邊界和波浪破碎時(shí)通常采用近似方法，引入過(guò)多可調(diào)參數(shù)，這不但給實(shí)際應(yīng)用帶來(lái)不便，某種程度上也加劇了模型的不穩(wěn)定。完全非線性淺水方程和Boussinesq水波方程同屬淺水方程范疇，忽略高階非線性和色散項(xiàng)后，Boussinesq水波方程可退化為完全非線性淺水方程。完全非線性淺水方程屬于典型雙曲守恒律方程，這類方程有許多基于高分辨率有限體積格式的高精度解法[10]，且多以黎曼間斷解問(wèn)題為理論基礎(chǔ)。而海岸水-陸動(dòng)邊界問(wèn)題本質(zhì)上屬于黎曼問(wèn)題，正適合采用完全非線性淺水方程進(jìn)行描述。此外，由于波浪破碎時(shí)Boussinesq方程中高階非線性項(xiàng)和色散項(xiàng)可以忽略不計(jì)[10]，破碎波仍可利用完全非線性淺水方程能夠捕捉間斷這一特點(diǎn)進(jìn)行計(jì)算，因此采用高分辨率有限體積格式求解Boussinesq類水波方程具有重要的意義。鑒于Boussinesq類水波方程的表達(dá)式較為復(fù)雜，目前以發(fā)展有限體積/有限差分混合格式為主，這類數(shù)值格式具備穩(wěn)定性強(qiáng)、處理波浪破碎和海岸水-陸動(dòng)邊界簡(jiǎn)單的特點(diǎn)[11-13]。

有限體積法的數(shù)值實(shí)現(xiàn)重點(diǎn)體現(xiàn)在網(wǎng)格界面處數(shù)值通量的計(jì)算，通常采用迎風(fēng)格式和中心格式。高階迎風(fēng)格式精度高，但需要精確或者近似求解黎曼問(wèn)題，構(gòu)造過(guò)程相對(duì)繁瑣，計(jì)算量大，對(duì)于實(shí)際應(yīng)用而言，較準(zhǔn)地求解黎曼解不是一件易事[10]；中心格式構(gòu)造簡(jiǎn)單，但容易出現(xiàn)數(shù)值震蕩。近幾年出現(xiàn)的中心迎風(fēng)格式兼具迎風(fēng)格式精確性和中心格式簡(jiǎn)單性優(yōu)點(diǎn)。中心迎風(fēng)格式無(wú)需求解復(fù)雜的黎曼問(wèn)題，程序?qū)崿F(xiàn)簡(jiǎn)單，目前采用中心迎風(fēng)格式求解Green-Naghdi水波方程中數(shù)值通量的研究并未見(jiàn)諸刊物，其精度和穩(wěn)定性等尚需進(jìn)一步考察。

文中將中心迎風(fēng)有限體積格式應(yīng)用于Green-Naghdi水波方程數(shù)值通量的求解，同時(shí)為了處理近岸波浪計(jì)算常遇到的水-陸動(dòng)邊界現(xiàn)象，采用靜壓重構(gòu)技術(shù)保證計(jì)算的穩(wěn)定性和水深h的非負(fù)性。該混合格式充分利用有限體積方法彌補(bǔ)有限差分方法在數(shù)值計(jì)算穩(wěn)定性和處理波浪破碎及水-陸動(dòng)邊界問(wèn)題的不足。通過(guò)模擬孤立波在常水深和斜坡海岸上的傳播及規(guī)則波跨越潛堤傳播實(shí)驗(yàn)對(duì)所建模型驗(yàn)證，數(shù)值計(jì)算結(jié)果與解析解和實(shí)驗(yàn)數(shù)據(jù)及其他模型計(jì)算所得結(jié)果進(jìn)行對(duì)比，結(jié)果表明，中心迎風(fēng)格式在程序編制難易程度、計(jì)算精度和計(jì)算效率等方面具有一定的優(yōu)勢(shì)。

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

1.1 控制方程

一維全非線性弱色散性Green-Naghdi方程表達(dá)式[6,15]：

其中，d=η+h為總水深，η為波面，h為靜水水深，u為水深平均速度，zb為水底高程，g為重力加速度。下標(biāo)t和x分別表示變量對(duì)時(shí)間t和空間x的偏導(dǎo)數(shù)。式中，Φ和Ψ分別為：

對(duì)于平底海床，水底高程zb為常數(shù)，即zb對(duì)空間變量x的導(dǎo)數(shù)為0，因此方程(1)可以改寫(xiě)為：

本文所建的一維Green-Naghdi水波方程模型將采用方程(1)作為控制方程。當(dāng)Φ=Ψ=0時(shí)，方程(1)和(3)分別與非平底海床和平底海床的淺水方程相一致。

為便于使用有限體積格式，將上述控制方程改寫(xiě)為守恒形式：

其中，U為變量矢量，F(xiàn)為通量矢量，定義為：

式中：

為利用具有和諧性質(zhì)的高精度有限體積格式處理海岸水-陸動(dòng)邊界問(wèn)題引入水位ζ(=zb+d)。式(4)中S為源項(xiàng)，方便起見(jiàn)，將其分成三個(gè)組成部分，分別為水底坡度項(xiàng)Sh、水底摩擦項(xiàng)Sf和色散項(xiàng)Sd，即：

其中，水底摩擦項(xiàng)由常用的二次律公式給出τ=-fu|u|，其中，f為水底摩擦系數(shù)，通常取值范圍是0.000 1～0.01。式(7)中的色散項(xiàng)由下式給出：

1.2 方程空間離散與時(shí)間積分

將式(4)在整個(gè)控制體單元格上積分，得到下式：

式中：Ω為單元區(qū)域；Γ為單元邊界；nx為外法向量。

采用矩形網(wǎng)格單元，將計(jì)算域在空間、時(shí)間上均勻劃分單元網(wǎng)格，做如下離散xi=iΔx(i=1,…,I)，tn=nΔt(n=1,…,N)，其中Δx，Δt分別為空間、時(shí)間步長(zhǎng)。波面η和速度u均定義在控制體中心。在有限體積[xi-1/2,xi+1/2]×[tn,tn+1]內(nèi)應(yīng)用方程(9)，采用對(duì)時(shí)間的差分近似為對(duì)時(shí)間的微分，并用差分代替導(dǎo)數(shù)，可得：

式中，Uin為n時(shí)刻解U在單元Ii∈[xi-1/2,xi+1/2]內(nèi)的空間平均值，F(xiàn)i+1/2為單元邊界xi+1/2處的數(shù)值通量，F(xiàn)i-1/2為單元邊界xi-1/2處的數(shù)值通量，Si為單元Ii對(duì)應(yīng)的數(shù)值源項(xiàng)。

基于Kurganov和Petrova提出的中心迎風(fēng)格式[16]，式(10)中的單元界面數(shù)值通量可以通過(guò)下式計(jì)算：

顯然，中心迎風(fēng)格式的有限體積法構(gòu)造相對(duì)簡(jiǎn)潔并且不需要求解復(fù)雜的黎曼問(wèn)題，因此中心迎風(fēng)格式的程序編制更為簡(jiǎn)單，二維情況下，這一優(yōu)勢(shì)將更加明顯。

上述中心迎風(fēng)格式僅為一階精度，為了提高空間精度，采用四階高精度狀態(tài)差值方法(MUSCL)[17]對(duì)界面左右變量進(jìn)行重構(gòu)，重構(gòu)后的界面變量再按照上述步驟進(jìn)行。模型的時(shí)間積分采用具有TVD性質(zhì)的三階龍格-庫(kù)塔方法[10]：

為保證計(jì)算收斂，時(shí)間計(jì)算步長(zhǎng)由Courant-Friedrichs-Lewy(CFL)穩(wěn)定性限制條件確定：

其中，ν為常數(shù)，對(duì)于本文計(jì)算，取ν=0.2。

1.3 邊界條件與波浪破碎的處理

模型中涉及的邊界條件包括固壁邊界、造波邊界和水-陸動(dòng)邊界。固壁邊界條件通過(guò)在計(jì)算域兩端分別設(shè)置3個(gè)虛擬網(wǎng)格實(shí)現(xiàn)，虛擬網(wǎng)格的變量通過(guò)固壁處法向速度為零，切向速度導(dǎo)數(shù)為零賦值。計(jì)算中涉及孤立波和規(guī)則波的生成，前者通過(guò)在計(jì)算域內(nèi)給定初始波面和流速作為入射波浪，后者采用在連續(xù)方程中增加源項(xiàng)的方法生成所需波浪，且根據(jù)需要在計(jì)算域兩端設(shè)置海綿層避免波浪二次反射的問(wèn)題[14]。水-陸動(dòng)邊界條件是波浪計(jì)算的難點(diǎn)，本文也是充分利用有限體積法求解完全非線性淺水方程的優(yōu)勢(shì)，將水-陸動(dòng)邊界處理為間斷。具體方法是通過(guò)在狀態(tài)重構(gòu)過(guò)程中(MUSCL)采用靜壓重構(gòu)技術(shù)實(shí)施[18]，該技術(shù)的特點(diǎn)就是將水-陸動(dòng)邊界附近的海床視為隨時(shí)間做微小變動(dòng)的，通過(guò)在計(jì)算過(guò)程中對(duì)海床的調(diào)整達(dá)到較準(zhǔn)確地捕捉水-陸動(dòng)邊界的目的。與常用的薄層水體法相比，該方法具有更高的精度，且保證了重構(gòu)水深非負(fù)和水-陸動(dòng)邊界處的和諧性，這對(duì)計(jì)算的穩(wěn)定性和精度非常重要。具體執(zhí)行步驟如下：

第一步：利用已知變量(p,ζ,d)采用MUSCL方法重構(gòu)網(wǎng)格界面左右狀態(tài)變量：

在此過(guò)程中，引入一極小水深作為判斷網(wǎng)格干、濕狀態(tài)的闕值(d<10-4m)，如某網(wǎng)格水深小于此值認(rèn)定為干網(wǎng)格，界面流速直接賦值為零。同理，可得到界面右側(cè)的速度值。需要注意的是，當(dāng)某一個(gè)網(wǎng)格的相鄰網(wǎng)格有干網(wǎng)格出現(xiàn)時(shí)，該網(wǎng)格以及干網(wǎng)格內(nèi)的變量分布認(rèn)為常數(shù)。

第二步：確定界面處水底高程值并進(jìn)行非負(fù)水深重構(gòu)，如下：

同理將上述過(guò)程用于右側(cè)界面，可得到右側(cè)界面狀態(tài)變量取值。

第三步：對(duì)水-陸邊界附近的網(wǎng)格進(jìn)行調(diào)整，即定義：

然后對(duì)zb進(jìn)行修正如下：

將通過(guò)式(16)-(20)進(jìn)行的變量重構(gòu)過(guò)程用于數(shù)值通量的計(jì)算。為了獲得和諧解，對(duì)于方程中的底坡源項(xiàng)要做如下相應(yīng)處理[10]：

處理波浪破碎問(wèn)題時(shí)采用文獻(xiàn)[10]給出的方法，即當(dāng)波高與水深之比達(dá)到0.8時(shí)，認(rèn)為波浪發(fā)生破碎，Green-Naghdi水波方程退化為完全非線性淺水方程，將破碎波浪處理為間斷。

2 數(shù)值驗(yàn)證

針對(duì)上文所建立的混合格式模型，本節(jié)將通過(guò)幾個(gè)典型算例的數(shù)值模擬分別對(duì)其進(jìn)行驗(yàn)證。其中，算例一為常水深孤立波的長(zhǎng)距離傳播問(wèn)題，用于驗(yàn)證模型的精度；算例二為破碎和非破碎孤立波在斜坡海岸上的傳播問(wèn)題，用于檢驗(yàn)?zāi)Ｐ筒蹲剿?陸動(dòng)邊界和模擬破碎波浪的能力；算例三為規(guī)則波跨越潛堤傳播，用于檢驗(yàn)?zāi)Ｐ蛯?duì)強(qiáng)非線性波浪模擬的能力。需要指出，所有計(jì)算均未使用人工光滑器。

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

孤立波是波浪色散性和非線性平衡制約的典型代表，其在常水深水槽中長(zhǎng)距離傳播的數(shù)值模擬是檢驗(yàn)數(shù)值模型的有力工具。而且，對(duì)于一維Green-Naghdi水波方程，孤立波存在精確解[19]，可以用來(lái)驗(yàn)證模型數(shù)值格式的精確性。

在本算例中，計(jì)算域長(zhǎng)度450 m，靜水水深h0為1.0 m，在計(jì)算域內(nèi)給出孤立波解析解作為初始條件，波高H=0.6 m，孤立波從計(jì)算域的左側(cè)向右側(cè)傳播，初始時(shí)刻波峰位于x0=50 m處，網(wǎng)格尺寸Δx=0.10 m，忽略底摩擦的影響，且計(jì)算域兩端采用固壁邊界條件。圖1(a)給出了采用中心迎風(fēng)格式模擬的孤立波在t=0 s, 20 s, 40 s, 60 s和80 s時(shí)的波面狀態(tài)。從圖中可見(jiàn)，經(jīng)過(guò)長(zhǎng)距離的傳播，孤立波的波形和波高值均保持穩(wěn)定，幾乎未發(fā)生變化，表明所建立的數(shù)值格式?jīng)]有引入偽色散和數(shù)值耗散。圖1(b)為在時(shí)間t=80 s時(shí)，模型計(jì)算結(jié)果與解析解的比較，可以看出兩者高度吻合，說(shuō)明本模型具有較好的精確性。

圖1 孤立波的長(zhǎng)距離傳播Fig. 1 The propagation of a solitary wave over a long distance

圖2 孤立波爬坡實(shí)驗(yàn)海床Fig. 2 The sketch of experimental setup of solitary wave runup

2.2 孤立波在斜坡海岸上的爬坡

本節(jié)將用所建立的Green-Naghdi模型模擬孤立波在斜坡上的爬高實(shí)驗(yàn)[20]，這些實(shí)驗(yàn)結(jié)果被廣泛的用于檢驗(yàn)波浪爬高模型，實(shí)驗(yàn)海床如圖2所示，是一個(gè)坡度為1∶19.85的斜坡海岸。

首先進(jìn)行非破碎孤立波的爬坡模擬，孤立波的波高與水深之比H/h0=0.018 5，靜水水深h0=0.39 m。計(jì)算域的長(zhǎng)度為90 m，網(wǎng)格長(zhǎng)度Δx= 0.05 m，底摩擦系數(shù)取cf= 0.001。對(duì)計(jì)算結(jié)果進(jìn)行以下無(wú)因次處理：

圖3中繪出了H/h0=0.018 5時(shí)，無(wú)量綱時(shí)間t*= 30, 40, 50, 60, 70時(shí)波形空間分布圖，x*=0對(duì)應(yīng)初始岸線。將數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，如圖3(a)和圖3(b)所示，數(shù)值計(jì)算結(jié)果稍超前于實(shí)驗(yàn)結(jié)果，圖3(c)和圖3(d)較好地復(fù)演了孤立波爬坡及回落過(guò)程，而圖3(e)中水面與斜坡交界處較明顯的差別，這可能是因?yàn)轲ば杂绊懙慕Y(jié)果[7]。對(duì)于該組物理模型實(shí)驗(yàn)，Li等[7]采用基于有限元方法的Green-Naghdi方程模型進(jìn)行過(guò)數(shù)值模擬，計(jì)算結(jié)果也一并在圖3中給出，可見(jiàn)其與本文計(jì)算結(jié)果基本一致。該算例說(shuō)明模型能夠有效地模擬非破碎孤立波在海岸上的傳播、爬坡和回落過(guò)程中波浪的運(yùn)動(dòng)狀態(tài)。

其次模擬破碎波孤立波在斜坡海床上的爬高。設(shè)置孤立波波高與水深之比H/h0=0.3，其他條件不變。如圖4，取無(wú)量綱時(shí)間t*= 15, 20, 25, 30時(shí)數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果以及Li等[7]的數(shù)值計(jì)算結(jié)果進(jìn)行對(duì)比，波面同樣采用無(wú)量綱形式。從圖4(b)中可以看到，數(shù)值計(jì)算結(jié)果由于孤立波的提前破碎，造成波面的計(jì)算結(jié)果小于實(shí)驗(yàn)結(jié)果，且波前近乎垂直，但這并未影響波浪爬坡和回落的演化過(guò)程。數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果吻合程度高，而偏小于Li等的數(shù)值計(jì)算結(jié)果，但各時(shí)刻的運(yùn)動(dòng)特征模擬均吻合較好。這說(shuō)明本文所建立的Green-Naghdi模型能充分描述破碎孤立波傳播、爬坡及回落過(guò)程中典型的運(yùn)動(dòng)特征，有效驗(yàn)證了模型的精度，同時(shí)說(shuō)明了本模型具有較好的處理水-陸動(dòng)邊界及波浪破碎的能力。

圖3 H/h0=0.018 5時(shí)孤立波波面數(shù)值計(jì)算結(jié)果和實(shí)驗(yàn)數(shù)據(jù)對(duì)比Fig. 3 Comparison of simulated and experimental surfaces of a solitary wave with H/h0=0.018 5

圖4 H/h0=0.3時(shí)孤立波波面數(shù)值計(jì)算結(jié)果和實(shí)驗(yàn)數(shù)據(jù)對(duì)比Fig. 4 Comparison of simulated and experimental surfaces of a solitary wave with H/h0=0.3

2.3 規(guī)則波在潛堤海床上的傳播

規(guī)則波跨越潛堤傳播是非常復(fù)雜的波浪演化過(guò)程，涉及波浪的變淺、高次諧波產(chǎn)生釋放等，是檢驗(yàn)?zāi)Ｐ蜕⑿院头蔷€性綜合性能的有力工具。采用Dingemans[21]的規(guī)則波跨越潛堤傳播的實(shí)驗(yàn)來(lái)對(duì)數(shù)值模型檢驗(yàn)，實(shí)驗(yàn)海床設(shè)置在圖5中給出，海床兩端設(shè)置1.5倍波長(zhǎng)海綿層吸收波浪，造波機(jī)在x=0 m處，平底處?kù)o水水深為0.4 m，梯形壩的坡腳位于x=6.0 m處，梯形壩頂端長(zhǎng)為2.0 m，距離靜水面0.1 m，梯形潛堤前坡坡度為1∶20，后坡坡度為1∶10，根據(jù)Dingemans[21]的實(shí)驗(yàn)，在計(jì)算域內(nèi)設(shè)置10個(gè)浪高儀，位置分別為x=2 m, 4 m, 10.5 m, 12.5 m, 13.5 m, 14.5 m, 15.7 m, 17.3 m, 19 m和21 m。入射波波高H=0.02 m，周期T=2.02 s，計(jì)算網(wǎng)格長(zhǎng)度為Δx=0.025 m。

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

圖6給出了10個(gè)浪高儀位置上數(shù)值模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比。如圖6所示，規(guī)則波跨越堤頂傳播過(guò)程中，波浪受到潛堤的影響，波形發(fā)生明顯變化，潛堤前以及堤頂位置處波面的數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合程度高，說(shuō)明模型能夠準(zhǔn)確模擬規(guī)則波的傳播及跨越潛堤時(shí)的波浪狀態(tài)。但是，在堤后浪高儀(例如x≥15.7 m)波面時(shí)間序列圖中可以看到，數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果在波幅和相位上均有較明顯地差別，這是因?yàn)橐?guī)則波跨越潛堤后的波浪中包含了具有強(qiáng)色散性的高頻自由波，而本文所使用的Green-Naghdi方程僅具有弱色散性，對(duì)強(qiáng)色散性波浪的模擬效果較差。使用高級(jí)別Green-Naghdi方程，例如GN-5方程，可得到更為精確的數(shù)值結(jié)果[22](本文采用方程為GN-1型方程)。

圖6 規(guī)則波跨越潛堤傳播各浪高儀處的波面時(shí)間序列Fig. 6 Time series of surface elevations for regular waves passing over a submerged bar

圖6中同時(shí)給出Li等[7]的對(duì)該物理模型實(shí)驗(yàn)的數(shù)值計(jì)算結(jié)果，可見(jiàn)其與本文計(jì)算結(jié)果基本一致，但存在一些細(xì)節(jié)差別。在圖6(a)-(f)測(cè)點(diǎn)處，本文模型計(jì)算結(jié)果與Li等的計(jì)算結(jié)果吻合程度較高，而在圖6(g)-(j)測(cè)點(diǎn)處Li等波面幅值計(jì)算結(jié)果均偏高于本模型計(jì)算結(jié)果，盡管如此，兩者差異較小，而且相位吻合程度高，這進(jìn)一步證明了數(shù)值模型的精度與正確性。

3 結(jié) 語(yǔ)

將中心迎風(fēng)格式應(yīng)用于建立求解一維全非線性Green-Naghdi水波方程的混合數(shù)值格式。對(duì)守恒形式控制方程中的通量項(xiàng)，采用有限體積方法離散和積分，應(yīng)用中心迎風(fēng)格式求解界面處數(shù)值通量，界面左右變量則通過(guò)四階狀態(tài)插值方法重構(gòu)，同時(shí)結(jié)合靜壓重構(gòu)技術(shù)保證計(jì)算域內(nèi)的和諧性和重構(gòu)水深非負(fù)性。方程中剩余項(xiàng)采用有限差分方法進(jìn)行求解。應(yīng)用具有TVD性質(zhì)的三階龍格-庫(kù)塔方法進(jìn)行時(shí)間積分。針對(duì)典型算例進(jìn)行了數(shù)值模擬，結(jié)果表明所建立的數(shù)值模式具有穩(wěn)定性強(qiáng)、間斷捕捉、處理海岸水-陸動(dòng)邊界及波浪破碎方便等優(yōu)點(diǎn)。應(yīng)用典型算例對(duì)所建模型進(jìn)行了驗(yàn)證，數(shù)值計(jì)算結(jié)果與解析解、實(shí)驗(yàn)數(shù)據(jù)及其他模型計(jì)算結(jié)果吻合良好，綜合考慮程序編制難易程度、計(jì)算精度和計(jì)算效率等方面，中心迎風(fēng)格式值得在Green-Naghdi水波方程的求解方面進(jìn)行推廣和應(yīng)用。

[1] CHOI D Y, WU C H, YOUNG C C. An efficient curvilinear non-hydrostatic model for simulating surface water waves [J]. International Journal for Numerical Methods in Fluids, 2011, 66(9): 1093-1115.

[2] KIRBY J T, WEI G, CHEN Q, et al. FUNWAVE 1.0: fully nonlinear Boussinesq wave model-Documentation and user's manual [R]. University of Delaware, 1998.

[3] LYNETT P J． Nearshore wave modeling with high-order Boussinesq-type equations [J]． Journal of Waterway, Port, Coastal and Ocean Engineering, 2006, 132(1): 348-357．

[4] 李孟國(guó), 王正林, 蔣德才. 近岸波浪傳播變形數(shù)學(xué)模型的研究與進(jìn)展 [J]. 海洋工程, 2002, 20(4): 43-57. (LI Mengguo, WANG Zhenglin, JIANG Decai. A review on the mathematical models of wave transformation in the nearshore region [J]. The Ocean Engineering, 2002, 20(4): 43-57. (in Chinese))

[5] 周俊陶, 林建國(guó), 謝志華. 高階Boussinesq類方程數(shù)值求解及試驗(yàn)驗(yàn)證 [J]. 海洋工程, 2007, 25(1): 88-92. (ZHOU Juntao, LIN Jianguo, XIE Zhihua. Numerical simulation for higher-order Boussinesq type equations and experimental verifications [J]. The Ocean Engineering, 2007, 25(1): 88-92. (in Chinese))

[6] GREEN A E, NAGHDI P M. A derivation of equations for wave propagation in water of variable depth [J]. Journal of Fluid Mechanics, 1976, 78(02): 237-246.

[7] LI M, GUYENNE P, LI F, et al. High order well-balanced CDG-FE methods for shallow water waves by a Green-Naghdi model [J]. Journal of Computational Physics, 2014, 257: 169-192.

[8] 趙彬彬, 段文洋. 全非線性深水波的Green-Naghdi理論研究 [J]. 哈爾濱工程大學(xué)學(xué)報(bào), 2009,8: 860-866. (ZHAO Binbin, DUAN Wenyang. Research on the Green-Naghdi theory for fully nonlinear deep water wave [J]. Journal of Harbin Engineering University, 2009,8: 860-866. (in Chinese))

[9] 朱良生, 洪廣文. 任意水深變化Boussinesq型方程非線性波數(shù)值計(jì)算 [J]. 海洋工程, 2000,18(2): 29-37. (ZHU Liangsheng, HONG Guangwen. Numerical calculation of no nlinear wave using Boussinesq equation in water of arbitrarily varying depths [J]. The Ocean Engineering, 2000,18(2): 29-37. (in Chinese))

[10] 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(2): 36-51.

[11] 岳良毅, 汪學(xué)鋒, 尤云祥. 擴(kuò)展型Boussinesq水波方程的一種改進(jìn)有限體積/有限差分混合算法[J]. 水動(dòng)力學(xué)研究與進(jìn)展：A輯, 2015,5:477-484. (YUE Liangyi, WANG Xuefeng, YOU Yunxiang. A modified hybrid finite-volume/finite-difference scheme for extended Boussinesq-type model [J]. Chinese Journal of Hydrodynamics, 2015,5: 477-484. (in Chinese))

[12] 房克照, 孫家文, 劉忠波, 等. Boussinesq水波方程新型數(shù)值解法 [J]. 海洋工程, 2014, 32(2): 97-103. (FANG Kezhao, SUN Jiawen, LIU Zhongbo, et al. Novel numerical solution to Boussinesq wave equations [J]. The Ocean Engineering, 2014, 32(2): 97-103. (in Chinese))

[13] ERDURAN K S, ILIC S, KUTIJA V. Hybrid finite-volume finite-difference scheme for the solution of Boussinesq equations [J]. International Journal for Numerical Methods in Fluids, 2005,49(11): 1213-1232.

[14] 房克照, 王磊, 劉忠波, 等. 基于MUSTA格式的全非線性Boussinesq波浪傳播數(shù)值模型 [J]. 力學(xué)學(xué)報(bào), 2014,5: 647-654. (FANG Kezhao, WANG Lei, LIU Zhongbo, et al. A fully nonlinear Boussinesq model for wave propagation based on MUSTA scheme [J]. Chinese Journal of Theoretical and Applied Mechanics, 2014,5: 647-654. (in Chinese))

[15] LE MéTAYER O, GAVRILYUK S, HANK S. A numerical scheme for the Green-Naghdi model [J]. Journal of Computational Physics, 2010, 229(6): 2034-2045.

[16] KURGANOV A, PETROVA G. A second-order well-balanced positivity preserving central-upwind scheme for the Saint-Venant system [J]. Communications in Mathematical Sciences, 2007, 5(1): 133-160.

[17] YAMAMOTO S, KANO S, DAIGUJI H. An efficient CFD approach for simulating unsteady hypersonic shock-shock interference flows [J]. Computers & Fluids, 1998, 27(5): 571-580.

[18] WANG Y, LIANG Q, KESSERWANI G, et al. A 2D shallow flow model for practical dam-break simulations [J]. Journal of Hydraulic Research, 2011, 49(3): 307-316.

[19] LANNES D, MARCHE F. A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations [J]. Journal of Computational Physics, 2015, 282: 238-268.

[20] SYNOLAKIS C E. The runup of long waves [D]. California Institute of Technology, 1986.

[21] DINGEMANS M W. Comparison of computations with Boussinesq-like models and laboratory measurements[R]. Deltares (WL), 1994.

[22] ZHAO B B, DUAN W Y, ERTEKIN R C. Application of higher-level GN theory to some wave transformation problems [J]. Coastal Engineering, 2014, 83: 177-189.

A numerical scheme for one-dimensional and fully nonlinear Green-Naghdi wave equations

JIAO Zifeng1, YIN Jing1, 2, FANG Kezhao1, SUN Jiawen1, 2

(1. The State Key Laboratory of Coastal and Offshore Engineering, Ocean Engineering Joint Research Center of DUT-UWA, Dalian University of Technology, Dalian 116024, China; 2. National Marine Environmental Monitoring Center, State Oceanic Administration, Dalian 116023, China)

O353.2

A

10.16483/j.issn.1005-9865.2016.06.004

1005-9865(2016)06-0030-08

2016-02-01

國(guó)家自然科學(xué)基金資助項(xiàng)目 (51579034)；中國(guó)科學(xué)院海洋環(huán)流與波動(dòng)重點(diǎn)實(shí)驗(yàn)室開(kāi)放基金課題資助(KLOCW1502);遼寧省教育廳重點(diǎn)實(shí)驗(yàn)室基礎(chǔ)研究項(xiàng)目資助(LZ2015013)

焦子峰(1989-)，男，山東日照人，碩士研究生，從事海岸動(dòng)力學(xué)方面研究。E-mail：zfjiao@mail.dlut.edu.cn