劉忠波， 劉澤， 王鋮龍， 王彥， 房克照

(1.大連海事大學(xué)交通運(yùn)輸工程學(xué)院， 大連 116026； 2.大連理工大學(xué)海岸和近海國家重點(diǎn)實(shí)驗(yàn)室， 大連 116024)

Boussinesq型水波方程是一類時(shí)域水波模型，此類模型最顯著的基本特征是控制方程中變量不對垂向坐標(biāo)z求導(dǎo)，等同于將三維波浪問題降至水平二維問題，同時(shí)Boussinesq模型也無需追蹤自由面問題，從而大大降低了計(jì)算時(shí)間，這促進(jìn)了其在海岸波浪水動(dòng)力模擬研究的廣泛應(yīng)用。

考慮到模型在實(shí)際工程應(yīng)用的可能性，Boussinesq型水波方程的變量對空間坐標(biāo)求導(dǎo)一般不超過3，關(guān)于這一點(diǎn)，很多理論和數(shù)值應(yīng)用文獻(xiàn)均有陳述[1-4]。吳昕煒等[5]利用Boussinesq波浪(Boussinesq wave，BW)模型對波浪在實(shí)際地形中的傳播變形規(guī)律進(jìn)行了數(shù)值模擬研究，并利用相關(guān)試驗(yàn)結(jié)果進(jìn)行了驗(yàn)證，該模型最高空間導(dǎo)數(shù)為3。最高空間導(dǎo)數(shù)為2的弱非線性Boussinesq型水波方程由Zhao等[6]給出，該方程具有Padé(2，2)的色散性能，一些學(xué)者研究泥沙動(dòng)力學(xué)問題前選擇了這個(gè)模型開展相關(guān)的波浪水動(dòng)力研究，但弱非線性是該模型的主要限制。Liu等[2]、Chazel等[7]給出最高空間導(dǎo)數(shù)為2的多層三維Boussinesq水波方程；Liu等[8]進(jìn)一步分析了Liu等[2]多層方程的理論性能，多層模型顯著提高了方程的線性和非線性性能。作為代價(jià)，層數(shù)越多，計(jì)算時(shí)間也會相應(yīng)增加。Liu等[8]模擬水槽聚焦波算例的結(jié)果表明，3層(11個(gè)方程)比2層(9個(gè)方程)計(jì)算時(shí)間增約為23%。因此，選擇單層模型也具有研究價(jià)值，盡管在1%色散誤差下原方程的適用最大因次水深kh=3.8，但精確的色散關(guān)系僅在kh<1.23左右(0.1%誤差)。與色散性能更為優(yōu)良的多層模型不同，單層模型在實(shí)際工程中被選擇應(yīng)用的概率更大，因此有限水深時(shí)的性能更應(yīng)予以重視。

除色散性外，淺化和非線性等性能是推導(dǎo)Boussinesq型水波方程必須檢測的重要指標(biāo)[2,8-10]。以變淺性能為例，劉忠波等[11]對水平二維Boussinesq方程[9,12-13]中的速度用變換速度取代，并通過分析變淺性能獲取了相關(guān)參數(shù)值。事實(shí)上，變換速度中引入水深變化項(xiàng)使得控制方程含有的獨(dú)立色散參數(shù)和獨(dú)立變淺參數(shù)的個(gè)數(shù)相等，這是方程具有更優(yōu)良變淺性能的關(guān)鍵因素。與水平二維方程不同，三維Boussinesq型方程含垂向速度，確保了方程具有更優(yōu)良的非線性性能[2,7-8,14]。

為了平衡效率和模型精度，現(xiàn)選擇Liu等[2]給出最高空間導(dǎo)數(shù)為2的單層三維Boussinesq水波方程為研究對象，對其理論性能的改進(jìn)方面進(jìn)行初探，并建立相應(yīng)的數(shù)值模型，通過相關(guān)試驗(yàn)來驗(yàn)證理論模型的改進(jìn)。

1 三維Boussinesq型方程的改進(jìn)

Liu等[2]在無旋、無黏假設(shè)下，首先沿垂向方向?qū)⑺w分為N層，并規(guī)定每層中點(diǎn)水深處的速度為特征速度，在每一層中，對該特征速度關(guān)于z做泰勒展開，并利用連續(xù)方程和無旋假定，則可消除掉變量對垂向坐標(biāo)z的導(dǎo)數(shù)，可以得到一組以N層特征速度表達(dá)的Boussinesq型水波方程。其次，對N個(gè)特征速度變量引入L變換算子得到N個(gè)變換速度(又稱為偽速度或計(jì)算速度)，引入變換速度的核心目的是提高方程的色散性，最終推導(dǎo)出具有更高性能的三組N層Boussinesq型水波方程。選取最高導(dǎo)數(shù)為2的單層方程作為研究對象，方程中的主要變量為波面位移η、自由表面處水平和垂向速度uη和wη、靜止水位處的速度u10和w10以及變換速度uE和wE。

在自由表面處，精確的運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)邊界條件為

(1)

(2)

U=uη+wη?η

(3)

式中：U為定義在自由表面處速度勢的水平梯度；g為重力加速度；?為水平梯度算子。

在海底，邊界滿足條件為

ub·?h+wb=0

(4)

式(4)中：ub和wb分別為海底水平速度和垂向速度；h為靜止水位下的水深。

自由表面處速度與靜止水位處速度滿足的關(guān)系表達(dá)式為

(5)

從水底到靜止水平面的速度場可表達(dá)為

(z-zα)?zα?·uE+(z-zα)?wE+

(6)

(z-zα)?zα·?wE-(z-zα)?·uE-

(7)

將z=0和z=-h分別代入式(6)和式(7)可得到靜止水位處和海底處的速度。在式(6)和式(7)中，引入的變換速度(uE,wE)與靜止的中間水深處速度(uα,wα)之間的L變換表達(dá)式為

(8)

(9)

與式(8)和式(9)給出的L變換表達(dá)式不同，引入的變換如下。

uα=uE+α1(αh)2?(?·uE)+β1α2h?h?·uE

(10)

wα=wE+α2(αh)2?·(?wE)+β2α2h?h·(?wE)

(11)

式中：α1、α2為色散參數(shù)；β1、β2為變淺參數(shù)。

由此，速度場表達(dá)式為

(z-zα)?zα?·uE+(z-zα)?wE+

[-2α2(z-zα)αh+(z-zα)2]?zα?·

(?wE)+β1α2h?h?·uE+

β2α2h(z-zα)?h?2wE

(12)

(z-zα)?zα·?wE-(z-zα)?·uE-

[-2α1(z-zα)αh+(z-zα)2]?zα·

?(?·uE)+β2α2h?h·(?wE)-

(z-zα)β1α2h?h·?2uE

(13)

類似于文獻(xiàn)[10]，在水底方程中引入變淺系數(shù)β3，最終水底運(yùn)動(dòng)學(xué)方程可寫為

(14)

式(1)、式(2)、式(5)、式(14)以及靜止水位處速度與變換速度的表達(dá)式[z=0代入式(12)和式(13)]，構(gòu)成了新三維Boussinesq型水波方程，當(dāng)α1=α2=1/6，β1=β2=0時(shí)，則可以轉(zhuǎn)化為原來的方程。

2 方程的基本性能

在立面二維和常水深情況下，通過傅里葉分析，則可分析出方程的相速度、和差頻、速度分布特征等，分析的詳細(xì)過程參見文獻(xiàn)[10]。

2.1 相速度

方程的色散關(guān)系表達(dá)式可寫為

(15)

式(15)中：ω為波浪圓頻率；k為波數(shù)。

當(dāng)表達(dá)式(15)確定的相速度與解析相速度在0

圖1 無因次相速度Fig.1 Non-dimensional wave phase celerity

2.2 速度分布

為了表征速度分布的精度，采用Madsen等[14]給出的公式考察本文方程與解析解之間的誤差，公式為

(16)

式(16)中：us(0) 和ws(0)分別為z=0時(shí)的Stokes波線性解。圖2給出了改進(jìn)值與原參數(shù)下的速度分布誤差的比較，整體來看，改進(jìn)參數(shù)后的水平速度和垂向速度誤差更小。

圖2 速度分布誤差Fig.2 The error of the velocity profiles

2.3 不規(guī)則波中的和差頻

一般來說，對于不含垂向速度的Boussinesq水波方程，差頻性能是核心短板，關(guān)于這一點(diǎn)，文獻(xiàn)[9,12,15]中均有表述。當(dāng)選取色散參數(shù)α1=0.170，α2=0.158(改進(jìn)值)和α1=α2=1/6(原參數(shù))情況下的，Boussinesq方程的和差頻解析結(jié)果如圖3所示，其中，k1和k2是兩個(gè)不同頻率波浪的波數(shù)。無因次和頻和差頻均以對角線為對稱線，對角線上面部分為和頻，下面部分是差頻，圖3僅給出一半的圖形。改進(jìn)后，1%誤差下和頻、差頻適用范圍均超過了kh=4。1%誤差下原參數(shù)方程的和頻、差頻適用范圍僅約等kh=1.7和kh=2.2。這表明，本文方程具有更優(yōu)良的和差頻特性。此外，無論是原方程還是本文方程，在非線性性能方面特別是差頻性能方面均優(yōu)于傳統(tǒng)型的Boussinesq水波方程，以Gobbi等[16]的方程為例進(jìn)行說明，該方程精確到4階保留了所有非線性項(xiàng)，色散適用水深約為kh=6，其1%誤差下差頻的適用范圍僅為kh=1左右。

圖3 無因次和差頻Fig.3 Non-dimensional amplitudes of sub-harmonics and super-harmonics

2.4 變淺性能

變淺性能代表著線性波浪由深水到淺水演化過程中的淺化特征，在0

圖4 方程的變淺梯度Fig.4 The shoaling gradient of the present model

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

對Boussinesq型水波方程可構(gòu)建相應(yīng)的垂向二維模型，數(shù)值模型中對時(shí)間差分采用3階Adams-Bashforth格式和4階Adams-Moulton格式，分別用于預(yù)報(bào)和校正階段。當(dāng)預(yù)報(bào)和校正的所有變量的誤差控制在一定范圍內(nèi)，則當(dāng)前計(jì)算步結(jié)束。波浪的產(chǎn)生采用Hsiao等[17]的域內(nèi)造波法，并在兩邊界上采用2倍波長的海綿邊界層進(jìn)行消波。

為對比改進(jìn)前后模型的差別，選擇了對Luth等[18]潛堤上規(guī)則波傳播變形的試驗(yàn)數(shù)值模擬。先簡要分析波浪在潛堤上的演化過程，波浪的傳播變形是一個(gè)復(fù)雜的物理過程，其演化過程中主要包括波-波相互作用、變淺與反變淺、相位鎖定與解鎖甚至波浪破碎等。為了精準(zhǔn)模擬演化過程，數(shù)值模型應(yīng)具有足夠精度的色散性、變淺性和非線性性能。

圖5 計(jì)算波面位移與文獻(xiàn)[14]試驗(yàn)結(jié)果的比較Fig.5 The comparisons of the computed surface elevations with reference[14]

數(shù)值模擬中，選取了周期T=2.02 s，波高為0.02 m的工況，時(shí)間步長采用0.01 s，空間步長采用0.03 m。圖5所示為改進(jìn)前與改進(jìn)后的計(jì)算波面位移與實(shí)驗(yàn)結(jié)果的比較。改進(jìn)后的本文方程和原方程的模擬結(jié)果與實(shí)驗(yàn)結(jié)果吻合程度均較高，但在浪高儀位置x=19.0 m和x=21.0 m處，本文方程比原方程在3次峰波形的準(zhǔn)確性上更高。具有Padé (2, 2)色散關(guān)系的各種Boussinesq數(shù)值模型無法準(zhǔn)確捕捉到最后3個(gè)位置點(diǎn)的波面位移。具有Padé (4, 4)色散關(guān)系的方程，若變淺性和非線性性能足夠，則可以較為精確捕捉到這3個(gè)位置點(diǎn)的波面，倘若變淺性能不佳，則也無法精準(zhǔn)模擬這一過程[19]。

4 結(jié)論

通過引入帶參數(shù)的變換速度，導(dǎo)出了一組最高空間導(dǎo)數(shù)為2的三維Boussinseq型方程。并進(jìn)行了理論分析和潛堤上波浪傳播變形模擬研究，主要得出以下結(jié)論。

(1)本文方程比原文獻(xiàn)在線性(相速度和速度分布)、和差頻等方面均有不同程度的改善，特別是相速度與解析解在0

(2)以實(shí)驗(yàn)結(jié)果為標(biāo)準(zhǔn)，潛堤后方高次諧波釋放，原方程對三次峰的捕捉不夠精確，這反映出原方程在色散性精度不足，但值得注意的是，本文方程模擬的主峰與實(shí)驗(yàn)結(jié)果吻合程度在最后兩個(gè)浪高儀位置與實(shí)驗(yàn)結(jié)果的吻合程度略次于原方程，關(guān)于這一點(diǎn)，更廣泛的數(shù)值模擬驗(yàn)證將會在下一步研究中進(jìn)行。