趙紅軍, 焦影霞, 孔俊

(1.河海大學(xué) 港口海岸與近海工程學(xué)院, 江蘇 南京 210098)

強(qiáng)非線性和色散性Boussinesq方程數(shù)值模型檢驗(yàn)

趙紅軍1, 焦影霞1, 孔俊1

(1.河海大學(xué) 港口海岸與近海工程學(xué)院, 江蘇 南京 210098)

采用同位網(wǎng)格有限差分法，建立了強(qiáng)非線性和色散性Boussinesq方程數(shù)值計(jì)算模型。以穩(wěn)恒波Fourier近似解給定入射波邊界條件，對(duì)均勻水深深水和淺水域不同非線性的行進(jìn)波、緩坡地形上深水至淺水域的淺水變形波、以及緩坡和陡坡地形上的波浪水槽實(shí)驗(yàn)進(jìn)行了數(shù)值計(jì)算，并將計(jì)算結(jié)果與解析解、解析數(shù)值解以及實(shí)驗(yàn)值進(jìn)行了較為詳細(xì)的比較，從而檢驗(yàn)了模型的色散性、非線性以及不同底坡下非線性波的淺水變形性能。

Boussinesq方程；非線性；色散性；淺水變形

1 引言

Boussinesq型方程是研究近岸水域波浪傳播變形的一個(gè)有力工具，雖然均勻水深水域下的Boussinesq方程[1]早在1872年就已給出，此類方程的實(shí)際應(yīng)用卻始于1960年代變水深水域經(jīng)典Boussinesq方程[2]的導(dǎo)得。因經(jīng)典Boussinesq方程僅具弱色散性和弱非線性，且未考慮水底摩阻和近岸波浪破碎的影響，所以自20世紀(jì)80年代始針對(duì)Boussinesq方程的研究主要集中在方程色散性和非線性的提高以及向近岸破波區(qū)的推廣。

在色散性能方面，Witting[3]首先引入了Padé近似和偽速度的概念，有力推動(dòng)了Boussinesq型方程色散性能的改進(jìn)。受Padé近似的啟發(fā)，Madsen等[4]在動(dòng)量方程中引入系數(shù)待定的三階導(dǎo)數(shù)項(xiàng)，得到了色散關(guān)系改進(jìn)型Boussinesq方程。Nwogu[5]推導(dǎo)了以任意層水平速度為特征速度的Boussinesq方程。任意層速度的引入不但提高了方程的色散性能，而且優(yōu)化了波動(dòng)水質(zhì)點(diǎn)速度的垂向分布，這在Gobbi及其合作者[6-7]以兩層水平速度勢(shì)函數(shù)重構(gòu)特征速度的理論研究中有著較好的體現(xiàn)。上述Boussinesq方程的改進(jìn)多以提高方程偏微分的階數(shù)為代價(jià)(如高至五階[7])，這為數(shù)值求解工作帶來了極大困難，為此，Lynett和Liu[8-9]引入了垂向分層的概念，他們以兩層或多層水平速度為特征速度表達(dá)的方程，在提高色散性能的同時(shí)還達(dá)到了偏微分降階的目的(降低至三階)。后來，Chazel等[10]進(jìn)一步將兩層Boussinesq方程偏微分的階數(shù)降低至兩階；Liu和Fang[11]的工作進(jìn)一步拓寬了兩層Boussinesq方程的水深適用范圍。

在非線性性能方面，早期的改進(jìn)方式主要是摒棄弱非線性假設(shè)，保留與色散項(xiàng)同階的所有非線性項(xiàng)，由此產(chǎn)生了所謂的“完全非線性Boussinesq方程”[12—13]；也有學(xué)者引入隨空間和時(shí)間變化的任意水深層概念[13—14]，以優(yōu)化該類方程的非線性性能。盡管早期的改進(jìn)工作在某些方面一定程度上提高了方程的非線性性能，但在水深較大時(shí)方程的非線性精度還遠(yuǎn)不及線性特征。Agnon等[15]在非線性改進(jìn)方面做出了富有成效的工作：他們基于Laplace方程和水底邊界條件，通過Padé近似和算子運(yùn)算技術(shù)，得到了靜水面水平速度和垂向速度之間的關(guān)系；以自由面速度和自由面為特征變量，表示自由面運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)邊界條件；通過速度勢(shì)級(jí)數(shù)展式建立自由面速度變量和靜水面速度變量之間的關(guān)系，從而給出了另一類形式的Boussinesq方程(稱之為Agnon型Boussinesq方程)。此后，Madsen等[16]通過算子運(yùn)算技術(shù)引入任意層偽速度，進(jìn)一步優(yōu)化提高了該類Boussinesq方程的色散性和水質(zhì)點(diǎn)速度的垂向分布；Wang等[17]通過任意層速度、Zhang等[18]通過任意層偽速度，分別給出了顯含定常背景水流影響的Agnon型Boussinesq方程。

盡管Agnon型Boussinesq方程具有良好的色散性能、變淺性能和非線性性能，但是繁雜的表達(dá)形式及其包含的高階導(dǎo)數(shù)項(xiàng)為方程的數(shù)值求解帶來了困難。Fuhrman和Bingham[19]構(gòu)建了空間七點(diǎn)差分、時(shí)間四階Runge-Kutta向前步進(jìn)的計(jì)算格式，對(duì)以任意層偽速度為特征速度的方程[16]進(jìn)行了數(shù)值求解；張洪生等[20]通過空間七點(diǎn)差分、時(shí)間五階Runge-Kutta-England向前步進(jìn)的格式，數(shù)值求解了以任意層速度為特征速度的方程[17]，并在文獻(xiàn)[21]中通過若干算例對(duì)含背景流影響的、以任意層速度[17]和以任意層偽速度[18]為特征速度的方程進(jìn)行了較為詳細(xì)的數(shù)值檢驗(yàn)。以往通過規(guī)則波算例對(duì)Agnon型Boussinesq方程的數(shù)值檢驗(yàn)[18,20-21]多以Stokes波或橢圓余弦波理論給定入射波邊界條件，然而任何給定入射波邊界與Boussinesq方程體系不符的情況均會(huì)導(dǎo)致計(jì)算過程產(chǎn)生高頻寄生波[22]而影響數(shù)值檢驗(yàn)效果，為此，本文以穩(wěn)恒波的Fourier近似解[23]提供入射波邊界，對(duì)以任意層速度為特征速度的Agnon型Boussinesq方程[17]開展了進(jìn)一步的數(shù)值檢驗(yàn)工作：首先，采用張洪生等[20]提出的同位網(wǎng)格有限差分格式，建立了以任意層速度為特征速度的Boussinesq方程數(shù)值計(jì)算模型；然后，通過均勻水深情況下深水和淺水域不同非線性水波傳播變形的計(jì)算，檢驗(yàn)了模型的色散性能和非線性性能；通過緩變和陡變斜坡地形上不同非線性水波的淺水變形計(jì)算，檢驗(yàn)了模型的淺化性能及其關(guān)于陡坡地形的適用性；最后，給出了結(jié)論。

2 Boussinesq方程數(shù)值模擬模型

2.1 控制方程與計(jì)算格式

基于Agnon等[15]的思路，Wang等[17]以任意層水平和垂向速度為特征速度，推導(dǎo)了緩變地形上含背景水流影響的具有強(qiáng)非線性和頻散性能的Boussinesq方程，無流情況下的空間一維方程體系可表達(dá)為：

(1)

(2)

(3)

(4)

(5)

(6)

式(1)～(6)構(gòu)成的偏微分方程組可采用同位網(wǎng)格有限差分法[20]進(jìn)行數(shù)值求解。因?yàn)榉匠探M包含了五階空間導(dǎo)數(shù)項(xiàng)，為確保低階導(dǎo)數(shù)項(xiàng)離散的誤差階數(shù)高于高階導(dǎo)數(shù)項(xiàng)，所以對(duì)方程組中的空間導(dǎo)數(shù)項(xiàng)統(tǒng)一采用七點(diǎn)有限差分格式進(jìn)行離散；同時(shí)，為保證時(shí)間導(dǎo)數(shù)項(xiàng)的計(jì)算誤差不混淆方程組中的五階空間導(dǎo)數(shù)，時(shí)間層上采用具有五階精度的變步長(zhǎng)Runge-Kutta-England格式向前步進(jìn)。此外，需要說明的是，以空間七點(diǎn)有限差分格式離散偏微分方程組(3)～(6)將形成一個(gè)稀疏系數(shù)矩陣方程組，為降低矩陣方程的求解難度，采用變量重新排序的方法[20]，將稀疏系數(shù)矩陣方程組變換為帶寬為27的帶狀矩陣方程組，并以列選主元的高斯消去法進(jìn)行矩陣方程求解。

2.2 邊界條件和數(shù)值濾波

2.2.1 入射波邊界條件

內(nèi)部造波源函數(shù)法[24—25]是Boussinesq方程數(shù)值模型入射邊界條件常用的處理方法，該法在計(jì)算域內(nèi)部指定位置處產(chǎn)生目標(biāo)波，造波的準(zhǔn)確性取決于造波源函數(shù)的合理解析。入射邊界條件亦可給定，如采用Stokes波理論或橢圓余弦波理論給定[18,20—21]。較造波源函數(shù)法，依據(jù)波浪理論解析解給定入射波邊界不僅可降低邊界誤差，而且可以考察Boussinesq方程與已有波浪理論的符合性，這是因?yàn)槿魏胃鶕?jù)解析解給定的邊界與Boussinesq方程體系不符的情況均會(huì)導(dǎo)致數(shù)值計(jì)算過程中產(chǎn)生高頻寄生波[22]。考慮到式(1)～(6)構(gòu)成的Boussinesq方程適用于深水至淺水的強(qiáng)非線性情況，所以本文采用穩(wěn)恒波Fourier近似理論[23]給定入射波邊界條件，這是因?yàn)檩^Stokes波和橢圓余弦波理論，F(xiàn)ourier近似理論關(guān)于深水至淺水波強(qiáng)非線性水波的解析具有更高的精度。此外，為了避免給定邊界處的初始不穩(wěn)定對(duì)數(shù)值計(jì)算的影響，我們對(duì)根據(jù)Fourier近似理論給定的入射波邊界值G乘以一隨時(shí)間變化的緩沖函數(shù)：

GR=Gtanh(t/NT)，

(7)

式中，N是常數(shù)，計(jì)算時(shí)取為5。

2.2.2 出流邊界條件

出流邊界條件的處理采用海綿層消波結(jié)合Sommerfeld輻射邊界條件的方式進(jìn)行。

(8)

式中，d為海綿層網(wǎng)格節(jié)點(diǎn)至邊界的距離；ds為海綿層寬度，當(dāng)取為1～2倍波長(zhǎng)的長(zhǎng)度時(shí)即可獲得較好的消波效果；α是海綿層緩沖參數(shù)，計(jì)算時(shí)取為3.5。

在一維情況下，Sommerfeld輻射邊界條件表達(dá)為：

(9)

2.2.3 數(shù)值濾波

在以同位網(wǎng)格有限差分格式求解偏微分方程組時(shí)，因中心差分格式的奇偶失聯(lián)，數(shù)值計(jì)算過程中會(huì)產(chǎn)生高頻的格點(diǎn)振蕩，數(shù)值實(shí)驗(yàn)顯示：這些高頻數(shù)值波動(dòng)的振幅有隨水深和非線性的增強(qiáng)而逐漸增強(qiáng)的趨勢(shì)。為消除這種數(shù)值高頻波對(duì)數(shù)值計(jì)算的影響，本文采用Bogey和Christophe[27]給出的具有低頻散性和低耗散性的九點(diǎn)選擇性濾波器進(jìn)行濾除，如下：

(10)

3 數(shù)值模型檢驗(yàn)

3.1 均勻水深水域波浪傳播的數(shù)值計(jì)算與檢驗(yàn)

設(shè)置長(zhǎng)度為20個(gè)波長(zhǎng)、絕對(duì)水深h=1.4 m的均勻水深水槽。水槽左端為入射波邊界，右端為開邊界，并在開邊界前設(shè)置長(zhǎng)度為100dx(dx為空間步長(zhǎng))的海綿層進(jìn)行消波。通過多種計(jì)算組合，對(duì)深水至淺水域的不同非線性水波的傳播變形進(jìn)行數(shù)值計(jì)算，具體的計(jì)算組合見表1，表中：H和T分別為波高和波周期，L是根據(jù)穩(wěn)恒波Fourier近似理論計(jì)算得到的波長(zhǎng)；算例包括深水波(h/L=1.0)和淺水波(h/L=0.05)兩種情況，對(duì)于深水波，非線性參數(shù)波陡(H/L)最大為0.08；對(duì)于淺水波，非線性參數(shù)波高水深比(H/h)最大為0.6。各種情況下的計(jì)算參數(shù)(包括空間步長(zhǎng)dx和時(shí)間步長(zhǎng)dt)如表1所列；為確保得到穩(wěn)定的計(jì)算結(jié)果，每種情況的計(jì)算總時(shí)間均取為50T。

圖1a～d和圖2a～d分別給出了深水域(h/L=1.0)和淺水域(h/L=0.05)中不同非線性水波在不同時(shí)刻(t=49.25T、49.50T、49.75T和50T)時(shí)無量綱波面(2η/H)的計(jì)算結(jié)果，圖中亦給出了穩(wěn)恒波Fourier近似解[23]的結(jié)果，為較為清晰地比較，圖件僅對(duì)14L～16L這一區(qū)域進(jìn)行了展示。結(jié)果顯示：數(shù)值計(jì)算的波形規(guī)則，波動(dòng)穩(wěn)定，隨著非線性的增強(qiáng)，波面不再關(guān)于靜水面對(duì)稱，波峰變得尖陡，波谷變得平坦。比較數(shù)值計(jì)算結(jié)果和Fourier近似解可知：對(duì)于深水和淺水域中的線性和弱非線性波(圖1a，b和圖2a，b)二者吻合良好，尤其是對(duì)于線性波的情況(圖1a和圖2a)，即使是對(duì)水深波長(zhǎng)比為1.0的深水波，二者也幾乎完全一致，說明模型具有優(yōu)良的色散性能。隨著非線性的增強(qiáng)，數(shù)值解與Fourier近似解的差異逐漸明顯：一方面表現(xiàn)為數(shù)值解的振幅較Fourier近似解逐漸減小，這應(yīng)該與數(shù)值計(jì)算的空間網(wǎng)格分辨率難以解析強(qiáng)非線性水波中的高頻約束波有關(guān)；另一方面表現(xiàn)為數(shù)值解的位相較Fourier近似解逐漸超前，說明對(duì)于強(qiáng)非線性的情況數(shù)值模型存在夸大波速的現(xiàn)象。圖3a，b分別給出了深水和淺水域中波速(c)和x=15L處波峰高度(ηc)的計(jì)算誤差(Error)隨入射波陡(H/L)和入射波波高水深比(H/d)的變化(誤差定義為Error=(f1-f2)/f2×100：f={c,ηc}；下標(biāo)1、2分別代表模型計(jì)算值和Fourier近似解；波速和波峰高度的計(jì)算值系通過上跨零點(diǎn)法據(jù)x≤15L范圍內(nèi)的數(shù)值波面統(tǒng)計(jì)得到)，誤差分析顯示：對(duì)于深水和淺水域中的強(qiáng)非線性水波，數(shù)值模型的波速誤差分別為0.21%和0.78%，在傳播15L后，波峰高度的誤差也僅為-3.2%和-3.7%。

表1 均勻水深水域波浪傳播變形實(shí)驗(yàn)算例

圖2 淺水域(h/L=0.05)、不同非線性(H/h=0.01～0.6)水波、不同時(shí)刻(t=49.25T(虛線)、49.50T(點(diǎn)劃線)、49.75T(雙點(diǎn)劃線)和50T(長(zhǎng)虛線))的無量綱自由面(2η/H)計(jì)算結(jié)果(虛線)及與Fourier近似解(t=50T，實(shí)線)的比較Fig.2 Computed non-dimensional free surface displacement (2η/H) at t=49.25T (dash), 49.50T (dash-dot), 49.75T(dash-dot-dot) and 50T (long-dash) from the simulations of wave propagation in constant shallow depth (d/L=0.05). Solutions of Fourier approxima-tion method at t=50T are drawn as the solid lines in a-d

圖3 深水域(a)和淺水域(b)中波速(◇)和波峰高度(○)的計(jì)算值較Fourier近似解的誤差Fig.3 Relative errors between the computed results and the Fourier approximated solutions in constant deep depth (a) and shallow depth (b)

3.2 緩坡地形上波浪傳播的數(shù)值計(jì)算：與解析解和解析數(shù)值解的比較

為了檢驗(yàn)?zāi)Ｐ偷臏\水變形性能，對(duì)不同非線性水波自深水至淺水域的淺水變形過程進(jìn)行了數(shù)值計(jì)算。圖4給出了數(shù)值實(shí)驗(yàn)的計(jì)算范圍及地形示意：數(shù)值水槽長(zhǎng)為120 m，在x=0 m處設(shè)置坡度為1∶100的緩變斜坡，斜坡上部接平臺(tái)；緩坡前平底段長(zhǎng)度設(shè)置為14 m，水深(h0)為1.0 m；為保證數(shù)值計(jì)算過程不致因水深變淺誘導(dǎo)的波浪破碎而發(fā)散，平臺(tái)水深(he)因入射波高的不同而有所變化(表2)。水槽右端設(shè)置為輻射邊界，并在輻射邊界前布置長(zhǎng)度為100dx的海綿層進(jìn)行消波；左端為給定的入射波邊界: 波周期T0為1.13 s，對(duì)應(yīng)的水深波長(zhǎng)比h0/L0=0.5(L0=gT2/2π是根據(jù)線性波理論計(jì)算的深水波長(zhǎng))；波高(H0)在0.000 02～0.08 m(表2算例a～d)，對(duì)應(yīng)的入射波陡H0/L0=0.000 01～0.04；其中算例a用于評(píng)價(jià)模型關(guān)于線性波淺水變形的計(jì)算效果，算例b～d用于考察模型關(guān)于有限振幅波淺水變形的計(jì)算效果。數(shù)值計(jì)算的空間步長(zhǎng)dx和時(shí)間步長(zhǎng)dt的取值分別如表2所列，各種情況下的計(jì)算總時(shí)間均為150T0，以最后5T0時(shí)間內(nèi)的波面計(jì)算結(jié)果通過上跨零點(diǎn)法計(jì)算波高。

圖4 緩坡地形上波浪傳播變形數(shù)值實(shí)驗(yàn)地形示意Fig.4 Bottom profile for wave propagation on mild slope topography

表2 緩坡地形上波浪傳播變形的實(shí)驗(yàn)算例

Tab.2 Cases for the computation of wave propagation on mild-slope topography

算例H0/mT0/she/mH0/L0dx/mT0/dta0.000021.130.040.000010.040150b0.021.130.040.010.024150c0.041.130.070.020.030150d0.081.130.120.040.030150

圖5 不同入射波陡情況下無量綱波面(2η/H0)的淺水變形過程Fig.5 Profiles of the non-dimensional free surface displacement (2η/H0) during the shoaling process under different incident wave steepness

圖6 不同入射波陡情況下無量綱波高(H/H0)的數(shù)值計(jì)算結(jié)果及其與解析解和解析數(shù)值解的比較Fig.6 Variations of the computational non-dimensional wave height (H/H0) and comparisons with those of linear theory and Fourier approximations during the shoaling process under different incident wave steepness

圖7 Ting和Kirby[29](a)、以及Tsai等[30](b)的物理模型實(shí)驗(yàn)斷面示意Fig.7 Experimental bottom profiles conducted by Ting and Kirby[29](a) and Tsai et al.[30](b)

圖5a～d分別給出了t=150T0時(shí)刻時(shí)不同入射波陡(H0/L0=0.000 01～0.04)情況下無量綱波面(2η/H0)的沿程變化，由此可見：規(guī)則波在緩坡地形上傳播時(shí)的波形規(guī)則，波態(tài)穩(wěn)定；對(duì)于微幅入射波(圖5a)的情況，數(shù)值波面關(guān)于靜水位對(duì)稱，且隨著水深變淺，行進(jìn)波的振幅呈先減小后增大的變化；對(duì)于有限振幅波的情況(圖5b，c)，隨著水深的減小，行進(jìn)波的非線性作用漸強(qiáng)，表現(xiàn)為數(shù)值波面逐漸抬高，峰谷關(guān)于靜水位不再對(duì)稱，尤以斜坡末端平臺(tái)起始位置處波面的波峰最為尖陡，波谷最為平坦，這是非線性作用生成高頻諧波的表現(xiàn)。在平臺(tái)初始位置處，算例b～d的局地波高水深分別為0.78、0.71和0.76，已接近或略小于孤立波一階近似理論確定的淺水極限波高水深比0.78；自此之后，平臺(tái)上的行進(jìn)波有隨傳播距離的增大而漸小的趨勢(shì)(圖5b～d)。

圖6a～d對(duì)不同入射波陡情況下無量綱波高(H/H0)的計(jì)算結(jié)果分別與解析解和解析數(shù)值解進(jìn)行了比較，其中，微幅入射波(H0/L0=0.000 01)的計(jì)算結(jié)果與基于線性波理論的淺水變形解析解比較(圖6a)；有限振幅入射波(H0/L0=0.01～0.04)的計(jì)算結(jié)果與基于穩(wěn)恒波傅里葉近似理論的解析數(shù)值解[28]比較(圖6b～d)。由此可知：在微幅入射波情況下(圖6a)，數(shù)值結(jié)果與線性理論解析解吻合甚好，雖然計(jì)算格式的數(shù)值耗散使計(jì)算值有隨傳播距離的增加(水深的減小)而逐漸小于解析解的趨勢(shì)，但二者的最大誤差(坡頂處)僅為-0.97%；在有限振幅入射波情況下(圖6b～d)，數(shù)值結(jié)果與解析數(shù)值解的差異亦隨水深的減小逐漸增大，尤以近坡頂處的差異為最大，并且這一差異有隨入射波陡的增加而逐漸增大的趨勢(shì)，如：在d/L0=0.1時(shí)，算例b、c和d淺水變形系數(shù)計(jì)算值較解析數(shù)值解的誤差分別為-0.83%、-1.17%和-2.26%；在d/L0分別為0.025、0.04和0.07時(shí)，算例b、c和d的計(jì)算值較解析數(shù)值解的誤差分別為-2.45%、-5.77%和-8.62%。

有限振幅入射波情況下淺水變形系數(shù)計(jì)算值較解析數(shù)值解偏小的原因可能在于如下3個(gè)方面：一是數(shù)值濾波，它可能會(huì)濾除非線性作用產(chǎn)生的高頻諧波而造成計(jì)算值偏??；其二源于微分方程差分近似產(chǎn)生的數(shù)值耗散，這一耗散與空間網(wǎng)格分辨率有關(guān)，一般情況下高的空間分辨率會(huì)產(chǎn)生小的數(shù)值耗散；其三可能是Boussinesq方程非線性淺化性能的限制所致，因?yàn)楫吘狗匠淌?3)中用于控制淺水變形性能的系數(shù)b1和b3是通過線性化的Boussinesq方程與線性波淺水變形解析解的比較確定的[17]，本節(jié)微幅波淺水變形計(jì)算結(jié)果與解析解相符合(圖6a)也輔證了線性理論情況下b1和b3系數(shù)的合理性。在數(shù)值濾波方面，本文采用了具有低頻散性和低耗散性的九點(diǎn)選擇性濾波器，根據(jù)Bogey和Bailly[27]的研究，此濾波器可有效濾除波長(zhǎng)小于4.7dx的波動(dòng)，這意味著對(duì)于非線性波的計(jì)算，如果在基頻波長(zhǎng)范圍內(nèi)設(shè)置47個(gè)空間網(wǎng)格步長(zhǎng)，那么這樣的空間分辨率可解析至10倍頻約束諧波。因?yàn)閿?shù)值耗散和數(shù)值濾波均與空間分辨率有關(guān)，所以為深入了解有限振幅情況下淺水變形計(jì)算值偏小的原因，我們對(duì)算例b～d分別在粗空間網(wǎng)格步長(zhǎng)dx=0.04 m下開展了進(jìn)一步的數(shù)值計(jì)算工作，并將相關(guān)計(jì)算結(jié)果示于圖6b～d中。比較不同空間分辨率下的結(jié)果可知：(1)在水深波長(zhǎng)比d/L0>0.08時(shí)(距入射波邊界為49L0)，不同空間分辨率下的結(jié)果基本重合，這說明數(shù)值耗散不是造成計(jì)算結(jié)果偏小的原因，因?yàn)槿绻菙?shù)值耗散所致，那么在長(zhǎng)距離傳播后不同空間網(wǎng)格步長(zhǎng)下的計(jì)算結(jié)果應(yīng)略存差異；(2)隨著水深波長(zhǎng)比的進(jìn)一步減小(d/L0<0.08)，不同空間分辨率下的計(jì)算差異逐漸明顯，粗網(wǎng)格下的數(shù)值結(jié)果明顯小于細(xì)網(wǎng)格，這是空間分辨率的不足造成非線性作用產(chǎn)生的高頻約束諧波因數(shù)值濾波器的耗散作用所致，同樣也是平臺(tái)上行進(jìn)波的振幅隨傳播距離逐漸減小的原因。因?yàn)樵赿/L0<0.08時(shí)不同入射波陡情況下的計(jì)算結(jié)果之于空間網(wǎng)格分辨率的敏感性甚小，所以有限振幅入射波淺水變形系數(shù)計(jì)算結(jié)果較解析數(shù)值解偏小的原因不是模型的數(shù)值耗散和數(shù)值濾波所引起，而應(yīng)該是控制方程在非線性淺化性能方面的局限所致。

3.3 斜坡地形上波浪傳播的數(shù)值計(jì)算：與物理模型實(shí)驗(yàn)的比較

以Ting和Kirby[29](記為TK)、以及Tsai等[30](記為Ts)進(jìn)行的波浪水槽實(shí)驗(yàn)對(duì)模型在緩變和陡變斜坡上的淺水變形過程進(jìn)行檢驗(yàn)，圖7給出了相應(yīng)的實(shí)驗(yàn)情況：其中，TK實(shí)驗(yàn)斷面是坡度為1/35的緩坡地形，Ts實(shí)驗(yàn)斷面為斜坡平臺(tái)地形，斜坡坡度包括1/10、1/5和1/3三種情況，平臺(tái)高度為0.8 m。表3給出了模型實(shí)驗(yàn)的水深(h0)、入射波浪條件(H0、T0)、斜坡坡度(cotβ)、破碎位置(xb為相對(duì)于坡腳的水平距離)、破碎類型以及波浪數(shù)值模型計(jì)算的空間步長(zhǎng)和時(shí)間步長(zhǎng)，表中，Li是根據(jù)穩(wěn)恒波Fourier近似理論計(jì)算得到的入射波波長(zhǎng)，需要說明的是：Ts實(shí)驗(yàn)并未對(duì)破碎位置和破碎類型進(jìn)行報(bào)道，表中給出的破碎位置是最大波高測(cè)波桿所在的位置(實(shí)際破碎位置可能產(chǎn)生于此位置前后)，破碎類型依據(jù)Galvin[31]的研究經(jīng)由Irribarren數(shù)的判定得到。此外，因?yàn)楝F(xiàn)階段的數(shù)值模型尚未包含任何破碎機(jī)制，所以為使計(jì)算過程不致因波浪破碎而發(fā)散，模型自破碎位置起至水槽末端設(shè)置100dx長(zhǎng)度的海綿層。

圖8a和8b分別給出了數(shù)值計(jì)算得到的TK兩組實(shí)驗(yàn)的波面(η)沿程變化(零點(diǎn)起自坡腳)，為方便比較，該圖亦將各測(cè)波桿記錄的峰谷實(shí)驗(yàn)值示于其中。結(jié)果顯示：數(shù)值計(jì)算得到的波形規(guī)則，波態(tài)穩(wěn)定；因斜坡上水深的減小，行進(jìn)波的淺水變形效應(yīng)顯著，峰谷關(guān)于靜水面表現(xiàn)出強(qiáng)烈的不對(duì)稱，尤以破碎位置處最為明顯；在經(jīng)過破碎位置后，行進(jìn)波的振幅因海綿層的作用而逐漸減弱。比較破碎前波峰和波谷的數(shù)值計(jì)算結(jié)果和實(shí)驗(yàn)值可知，二者吻合良好，數(shù)值結(jié)果較好地反映了有限水深和淺水域中非線性水波淺水變形過程中的峰谷變化。

表3 斜坡地形上波浪傳播變形的物理模型實(shí)驗(yàn)算例和數(shù)學(xué)模型計(jì)算參數(shù)

圖8 波面(η)沿程變化計(jì)算值(實(shí)線)與實(shí)驗(yàn)值(實(shí)心圓點(diǎn))的比較Fig.8 Variations of the computed free surface displacement (η) and comparisons with the experimental data

圖9 波高(H)沿程變化的計(jì)算值(實(shí)線)與實(shí)驗(yàn)值(實(shí)心圓點(diǎn))的比較Fig.9 Variations of the computed wave height (H) and comparisons with the experimental results

圖9a～e分別給出了數(shù)值計(jì)算得到的TK和Ts 5組實(shí)驗(yàn)情況下波高的沿程變化，結(jié)果顯示：對(duì)于緩坡(1/35)TK1和TK2實(shí)驗(yàn)，斜坡上的波高隨水深的減小而緩慢增大；對(duì)于陡坡(1/10、1/5和1/3)Ts1、Ts2和Ts3實(shí)驗(yàn)，斜坡上的波高因水深的劇減而快速增強(qiáng)，斜坡前有一定程度的反射波生成，且反射波的振幅有隨坡度的增大而逐漸增強(qiáng)的趨勢(shì)。從最大波高的計(jì)算結(jié)果與實(shí)驗(yàn)值的比較可以看出：對(duì)于崩破波實(shí)驗(yàn)TK1，最大波高的發(fā)生位置和強(qiáng)度均與實(shí)驗(yàn)值吻合良好，且結(jié)合圖8a可知，此時(shí)波峰和波谷值亦與實(shí)驗(yàn)值相吻合；對(duì)于卷破波TK2以及Ts1～Ts3實(shí)驗(yàn)，數(shù)值計(jì)算的最大波高值均略低于實(shí)驗(yàn)值，陡坡Ts1～Ts3實(shí)驗(yàn)中的最大波高位置與實(shí)驗(yàn)值相吻合，緩坡TK2實(shí)驗(yàn)中的最大波高位置較實(shí)驗(yàn)結(jié)果略為超前，這可能與Boussinesq方程無法描述卷波破碎時(shí)的波面翻卷現(xiàn)象有關(guān)。

4 結(jié)論

采用同位網(wǎng)格有限差分法，建立了強(qiáng)非線性和色散性Boussinesq方程數(shù)值計(jì)算模型。以穩(wěn)恒波Fourier近似解提供入射波邊界條件，對(duì)均勻水深深水和淺水域的行進(jìn)波、以及緩變和陡變斜坡地形上不同非線性水波的淺水變形過程進(jìn)行了數(shù)值計(jì)算，數(shù)值結(jié)果與解析解、解析數(shù)值解和波浪水槽實(shí)驗(yàn)值的比較表明：(1)模型具有良好的色散性能和非線性性能，至少可對(duì)深水(水深波長(zhǎng)比1.0)強(qiáng)非線性(波高波長(zhǎng)比0.08)以及淺水(水深波長(zhǎng)比0.05)強(qiáng)非線性(波高水深比0.6)的行進(jìn)波進(jìn)行有效計(jì)算；(2)模型具有較好的淺水變形性能，可描述深水至淺水域弱非線性入射波的淺水變形過程；(3)盡管控制方程受限于緩變水深這一假設(shè)，但模型可對(duì)乃至1/3陡變斜坡上的淺水變形波進(jìn)行較為合理地預(yù)測(cè)。

[1] Boussinesq J. Theory of wave and swells propagated in long horizontal rectangular canal and imparting to the liquid contained in this canal[J]. Journal de Mathematiques Pures et Appliquees, 1872, 17(2): 55-108.

[2] Peregrine D H. Long waves on a beach[J]. Journal of Fluid Mechanics, 1967, 27: 815-827.

[3] Witting J M. A unified Model for the evolution of nonlinear water waves[J]. Journal of Computational Physics, 1984, 56: 203-236.

[4] Madsen P A, Murray R, S?rensen O R. A new form of the Boussinesq equations with improved linear dispersion characteristics[J]. Coastal Engineering, 1991, 15(4): 371-388.

[5] Nwogu O. Alternative form of Boussinesq equations for nearshore wave propagation[J]. Journal of Waterway, Port, Coastal, and Ocean Engineering, 1993, 119(6): 618-638.

[6] Gobbi M F, Kirby J T. Wave evolution over submerged sills: test of a high-order Boussinesq model[J]. Coastal Engineering, 1999, 37(1): 57-96.

[7] Gobbi M F, Kirby J T, Wei G. A fully nonlinear Boussinesq model for surface waves.Part 2. Extension to O(kh)4[J]. Journal of Fluid Mechanics, 2000, 405: 181-210.

[8] Lynett P, Liu P L-F. Linear analysis of the multi-layer model[J]. Coastal Engineering, 2004, 51: 439-454.

[9] Lynett P, Liu P L-F. A two-layer approach to wave modeling[J]. Proceeding of the Rval Society A, 2004, 460: 2637-2669.

[10] Chazel F, Benoit M, Ern A, et al. A double-layer Boussinesq-type model for highly nonlinear and dispersive waves[J]. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2009, 465: 2319-2346.

[11] Liu Z B, Fang K Z. Two-layer Boussinesq models for coastal water waves[J]. Wave Motion, 2015, 57: 88-111.

[12] Wei G, Kirby J T, Grilli S T, et al. A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[J]. Journal of Fluid Mechanics, 1995, 294: 71-92.

[13] Gobbi M F, Kirby J T, Wei G. A fully nonlinear Boussinesq model for surface waves, Part 2: Extension to O(kh)4[J]. Journal of Fluid Mechanics, 2000, 405: 18l-210.

[14] Zou Z L. A new form of higher order Boussinesq equations[J]. Ocean Engineering, 2000, 27(5): 557-575.

[15] Agnon Y, Madsen P A, Sch?ffer H A. A new approach to high-order Boussinesq model[J]. Journal of Fluid Mechanics, 1999, 399: 319-333.

[16] Madsen P A, Bingham H B, Liu H. A new Boussinesq method for fully nonlinear waves from shallow to deep water[J]. Journal of Fluid Mechanics, 2002, 462: 1-30.

[17] Wang Yalin, Zhang Hongsheng, Miao Guoping, et a1. A new approach to high order Boussinesq type equations with ambient currents[J]. China Ocean Engineering, 2005, 19(1): 49-60.

[18] Zhang Hongsheng, Wang Weiyuan, Feng Wenjing, et a1. A numerical model for nonlinear wave propagation on non-uniform current[J]. China Ocean Engineering, 2010, 24(1): 15-28.

[19] Fuhrman D R, Bingham H B. Numerical solutions of fully non-linear and highly dispersive Boussinesq equations in two horizontal dimensions[J]. International Journal of Numerical Methods in Fluids, 2004, 44(3): 231-225.

[20]張洪生, 馮文靜, 王亞玲,等. 非線性波傳播的新型數(shù)值模擬模型及其實(shí)驗(yàn)驗(yàn)證[J]. 海洋學(xué)報(bào), 2007, 29(4): 137-147.

Zhang Hongsheng, Feng Wenjing, Wang Yaling, et al. A new approach to numerical simulation of nonlinear wave propagation and its experimental verification[J]. Haiyang Xuebao, 2007, 29(4): 137-147.

[21] Zhang Hongsheng, Wang Weijuan, Feng Wenjing, et al. Tests and applications of a Boussinesq model with ambient current[J]. Journal of Hydrodyanmics, 2010, 22(4): 526-536.

[22] McCabe M. Modelling nearshore waves, runup and overtopping[D]. Manchester: University of Manchester, 2011.

[23] Zhao Hongjun, Song Zhiyao, Li Ling, et al. On the Fourier approximation method for steady water waves[J]. Acta Oceanologca Sinica, 2014, 33(5): 37-47.

[24] Larsen J, Dancy H. Open boundaries in short wave simulations: a new approach[J]. Coastal Engineering, 1983, 7: 285-297.

[25] Schaffer H A, Sorensen O R. On the internal wave generation in Boussinesq and mild-slope equations[J]. Coastal Engineering, 2006, 53(4): 319-323.

[26] Fuhrman D R, Madsen P A, Bingham H B. Numerical simulation of lowest-order short-crested wave instabilities[J]. Journal of Fluid Mechanics, 2006, 563: 415-441.

[27] Bogey C, Christophe B. A family of low dispersive and low dissipative explicit schemes for flow and noise computations[J]. Journal of Computational Physics, 2004, 194(1): 194-214.

[28] Zhao H J, Kong J, Song Z Y. Shoaling of nonlinear surface water waves on depth uniform current[C]//Processing 35th Congress International Assoc Hydraul Research. Chengdu, 2013.

[29] Ting F C K, Kirby J T. Observation of undertow and turbulence in a laboratory surf zone[J]. Coastal Engineering, 1994, 24(1/2): 51-80.

[30] Tsai C P, Chen H B, Hwung H H, et al. Examination of empirical formulas for wave shoaling and breaking on steep slopes[J]. Ocean Engineering, 2005, 32(3/4): 469-483.

[31] Galvin C J. Breaker type classification on three laboratory beaches[J]. Journal of Geophysical Research, 1968, 73(12): 3651- 3659.

Numerical validation of a Boussinesq-type model for highly nonlinear and dispersive waves

Zhao Hongjun1, Jiao Yingxia1, Kong Jun1

(1.CollegeofHarbor,CoastalandOffshoreEngineering,HohaiUniversity,Nanjing210098,China)

In the present work a numerical highly nonlinear and dispersive Boussinesq-type model is developed based on a non-staggered finite difference technique. With the Fourier approximation method providing the incident wave boundary condition, the model is applied and verified against a set of three test cases for which analytical, numerical or experimental reference results are available: (1) propagation of linear and nonlinear periodic waves on deep and shallow depth, (2) shoaling of linear and nonlinear regular waves from deep to shallow water on a mild slope, and (3) transformation of regular waves on a mild slope and on steep slopes. Comparisons of the numerical results with the analytical, numerical and experimental ones confirm the capabilities of the model for the predictions of highly nonlinear and dispersive waves and for the computations of nonlinear wave shoaling on different slopes.

Boussinesq equation; nonlinearity; dispersion; shoaling

10.3969/j.issn.0253-4193.2017.05.002

2016-10-14；

2017-01-15。

江蘇省自然科學(xué)基金青年基金項(xiàng)目 (BK20130827)；交通部重點(diǎn)科技項(xiàng)目(2015328521280)；水利部公益性科研專項(xiàng)(201501010)。

趙紅軍(1980—)，男，天津市薊縣人，博士，副教授，主要從事水波動(dòng)力學(xué)理論與應(yīng)用研究。E-mail:loyhg@hhu.edu.cn

P731.22

A

0253-4193(2017)05-0010-12

趙紅軍, 焦影霞, 孔俊. 強(qiáng)非線性和色散性Boussinesq方程數(shù)值模型檢驗(yàn)[J]. 海洋學(xué)報(bào), 2017, 39(5): 10-21,

Zhao Hongjun, Jiao Yingxia, Kong Jun. Numerical validation of a Boussinesq-type model for highly nonlinear and dispersive waves [J]. Haiyang Xuebao, 2017, 39(5): 10-21, doi：10.3969/j.issn.0253-4193.2017.05.002