三維Boussinesq型水波方程的改進

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

(1.大連海事大學交通運輸工程學院， 大連 116026； 2.大連理工大學海岸和近海國家重點實驗室， 大連 116024)

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

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

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

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

1 三維Boussinesq型方程的改進

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

在自由表面處，精確的運動學和動力學邊界條件為

(1)

(2)

U=uη+wη?η

(3)

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

在海底，邊界滿足條件為

ub·?h+wb=0

(4)

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

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

(5)

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

(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變換表達式為

(8)

(9)

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

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ù)。

由此，速度場表達式為

(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)

類似于文獻[10]，在水底方程中引入變淺系數(shù)β3，最終水底運動學方程可寫為

(14)

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

2 方程的基本性能

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

2.1 相速度

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

(15)

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

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

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

2.2 速度分布

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

(16)

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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