雙層Boussinesq水波方程速度公式的修正

劉忠波 ， 韓青亮， 任雙雙， 王 彥， 房克照

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

波浪是近岸水域重要的水動力之一，波浪力是海洋(海岸)工程結(jié)構(gòu)設計中必須考慮的環(huán)境荷載，精準預報這一波浪條件是工程設計的前提.水平波浪力是波浪壓力沿水深方向的積分，與波浪波面下的速度分布息息相關(guān)，因此對速度場的研究具有重要的學術(shù)和工程應用價值.

波浪速度的研究多采用理論分析、數(shù)值模型模擬和物理模型試驗模擬等方法，其中數(shù)值模型模擬省時且成本低，是最常用的方式.作為一種典型的勢流模型，多數(shù)Boussinesq方程將復雜的三維水波問題簡化為二維問題，大大降低了模型求解難度，促使這類方程在海岸波浪水動力研究中得以長足發(fā)展.Boussinesq模型的理論發(fā)展經(jīng)歷了從弱非線性到強非線性等過程，最新方程的適用水深已得到大幅度拓展[1-2]，甚至在Liu等[2]的研究中已完全擺脫了水深限制.文獻[1]首次對兩層水體中間位置的水平和垂直方向速度沿垂向坐標z做泰勒(Taylor)展開，進而利用計算速度取代中間位置處的速度，最終推導出最大空間導數(shù)為3的雙層Boussinesq方程，該方程在1%誤差內(nèi)的最大適用水深達kh=53.1(k為波數(shù)，h為靜水深)；文獻[2]進一步將文獻[1]拓展成多層，并給出空間導數(shù)為2和5的多層Boussinesq水波方程，4層空間導數(shù)為5的方程在1%誤差下的適用水深達kh=7 600.關(guān)于Boussinesq水波方程在理論性能、數(shù)值格式和數(shù)值應用的研究進展可參見文獻[3-5].Madsen等[6]的最新研究認為多數(shù)高階Boussinesq水波方程因色散精度不足，存在波谷不穩(wěn)定的缺點，并證明了文獻[2]最高導數(shù)為3的4層水波方程的穩(wěn)定范圍最大.

Liu等[7]分析了最高空間導數(shù)為2的多層Boussinesq水波方程的線性和非線性性能，1%誤差內(nèi)，2層模型色散適用水深(相速度)可達kh=19.7，而沿水深分布的水平速度的適用水深僅為kh=5.1，前者是后者的3.86倍；3層和4層色散適用水深分別是其水平速度的3.65倍和2.57倍，可見多層模型速度的適用水深遠小于相速度，這一現(xiàn)象也存在于多數(shù)常用Boussinesq水波方程中.林鵬程等[8]針對不同非線性、不同水深情況和線性入射條件，研究了最高空間導數(shù)為3的單層Boussinesq水波方程垂向分布的速度特征，得出的速度分布與Stokes線性波、二階波和三階波解析解存在一定的差異；劉必勁等[9]采用單層Boussinesq水波方程模擬聚焦波和穩(wěn)態(tài)波，發(fā)現(xiàn)單層Boussinesq水波方程能夠勝任波面的模擬，但無法給出精確的速度輪廓.文獻[8-9]采用的單層Boussinesq水波方程是文獻[1]雙層模型的簡化版，相關(guān)模擬結(jié)果也說明單層Boussinesq水波方程的速度精度遠小于色散精度.

最高空間導數(shù)為2的N層Boussinesq水波方程[2]的色散關(guān)系是Stokes波色散關(guān)系式的Padé[4N-2, 4N]型展開，對色散關(guān)系做泰勒展開則可確定相速度精確到kh的8N-4～8N次冪，而每一層速度沿水深分布的表達式中z的冪次僅為2，因此N層方程的速度最多可精確到kh的4N次冪，可見色散精度遠高于速度精度.保持方程表達形式不變的情況下，如何彌補這一結(jié)果帶來的差異，歷史文獻沒有給出明確答案.因此，本文選擇文獻[2]中最高導數(shù)為2的雙層Boussinesq水波方程作為研究對象，嘗試提高速度表達式中z的冪次，以獲得精度更高的速度場.

1 速度公式的修正

文獻[2]將垂向水體分為2層，推導了最高導數(shù)為2的雙層Boussinesq方程.在立面二維情況下，速度表達式分為從自由表面到靜水位、從靜水位到連接面以及連接面到水底共3段.

從靜水位到連接面的速度場為

(1)

(2)

從連接面到水底的速度場為

(3)

(4)

從自由表面到靜水位的速度場為

(5)

(6)

式中：u10,w10分別為靜水位處的水平和垂向速度.

利用式(1)～(6)可以計算整個水體內(nèi)的速度，而控制方程的表達形式保持不變[2,7]，以下簡稱為文獻[2]公式.為提高速度精度，在表達式(1)～(4)中引入乘以常數(shù)的3階項，速度表達式可改寫為

(7)

(8)

(9)

(10)

式中：β為常系數(shù)，其取值需通過與速度解析解做優(yōu)化獲取.式(7)～(10)中最后一項是3階導數(shù)項的修正，此3階項可在方程推導過程中直接得到.當速度位于z=zα1，z=zα2時，該項自動為0；當不在這兩處時，乘以β后相當于修正.式(5)～(10)簡稱為改進公式.

Sun等[10]研究發(fā)現(xiàn)，當色散參數(shù)α1取值范圍在0.13～0.25時，方程的相速度、變淺性、二階非線性和三階非線性等性能在0

(Fw(kh))2]d(kh)

(11)

水平速度誤差和垂直速度誤差采用Liu等[2]給出的表達式

(12)

式中：us(0)和ws(0) 為靜水位處的水平和垂向速度；us(z)和ws(z)為Stokes線性波速度場.取 0

將本文的改進公式與文獻[2]中的公式計算結(jié)果比較，如圖1和圖2所示. 由圖可見，誤差在1%內(nèi)，改進公式的水平速度和垂向速度的最大適用水深分別為kh=7.34和kh=7.83，均超過文獻[2]中公式的最大適用水深kh=5.1(水平速度)和kh=4.5(垂向速度)，表明改進公式在理論層面上是有效的.

圖1 水平速度的積分誤差Fig.1 Integrated error of horizontal velocity along water depth

圖2 垂向速度的積分誤差Fig.2 Integrated error of vertical velocity along water depth

為進一步對比文獻[2]中公式和改進公式的差異，給出kh=3.0，6.0，9.0時，沿垂向的速度分布與Stokes線性波速度解析解，如圖3所示.伴隨水深kh的增大，改進公式與文獻[2]公式的差異也會越來越明顯，解析解與本文改進公式的結(jié)果吻合程度更高.

圖3 水平和垂向速度與Stokes線性波解析解比較Fig.3 Comparisons of velocity profiles between numerical results and Stokes linear wave solution

2 數(shù)值模型及驗證

2.1 數(shù)值模型

本文數(shù)值模型依賴的方程和文獻[2]中原始方程相同，對應的數(shù)值模型在時間差分格式上采用了混合4階Adams-Bashforth-Moulton的時間步進格式，空間差分格式采用中心差分格式，邊界點上采用偏差公式.根據(jù)實際計算情況，可考慮在入射邊界上采用線性邊界造波技術(shù)、邊界松弛造波技術(shù)或在域內(nèi)采用內(nèi)部造波技術(shù)，在傳出邊界上采用海綿邊界層吸收波浪.具體過程可參考文獻[2,7,11].

2.2 與流函數(shù)波浪速度場的比較

圖3僅給出了線性波條件的速度解析對比，對于較強非線性情況下的速度場是否能夠勝任，還需通過數(shù)值模擬加以印證.以流函數(shù)波浪速度為比較對象，利用解析解驗證改進前后計算結(jié)果的差異，進一步驗證本文公式的有效性.

設計靜水深h分別為50 m和70 m的兩個工況，流函數(shù)波浪周期為T=6s、波高H=5 m.數(shù)值模擬中，L為流函數(shù)波浪對應的非線性波長，計算區(qū)域采用10L，空間步長采用L/32，時間步長采用0.05 s.改進公式和文獻[2]公式計算結(jié)果與解析解的比較如圖4所示，其中c為流函數(shù)波浪的相速度.與文獻[2]公式相比，改進公式計算得到的速度場與解析解吻合程度更佳，特別是在水深為70 m時，波浪流函數(shù)的波長為60.171 m，反算得出kh=7.31.結(jié)合圖1，當kh=7.31時，文獻[2]公式的水平速度積分誤差為1.83%，遠大于改進公式的積分誤差0.99%.數(shù)值結(jié)果表明，改進的速度場計算結(jié)果明顯優(yōu)于文獻[2]公式的結(jié)果.

圖4 計算波峰面下的速度場與波浪流函數(shù)解析速度場比較Fig.4 Comparison of computed velocity profile with analytical solution of wave stream function under wave crest

2.3 聚焦波最大波峰面下的水平速度剖面

Baldock等[12]進行了深水聚焦波演化試驗的研究，將最大和最小周期分為29份，每個頻率的波幅取值相同，其中B組為寬譜，T=0.6～1.4 s，kh=1.568～7.825；D組為窄譜，T=0.8～1.2 s，kh=2.026～4.403.

為了模擬聚焦波生成的全過程，在入射邊界處引入與試驗類似的線性疊加方式，將整個數(shù)值計算域設為24 m，空間步長為0.04 m，時間步長為 0.01 s，末端設置4 m的海綿邊界層以消波，總模擬時間為40 s，設定聚焦波峰面為55 mm的B組與D組進行數(shù)值模擬.Liu等[7]已建立了2層、3層和4層水波模型，模擬了文獻[12]中D55和B55這兩組非線性較強的工況，因波面計算結(jié)果仍與文獻結(jié)果保持一致，故不再重復給出相關(guān)結(jié)果，本文重點比較改進公式和文獻[2]公式的計算結(jié)果與這兩組工況的試驗結(jié)果，如圖5～6所示.

由圖可見，改進公式優(yōu)于文獻[2]公式的計算結(jié)果.D55窄譜情況下，波數(shù)范圍在kh= 2.026～4.403之間，結(jié)合圖1結(jié)果可知，文獻[2]中水平速度在1%誤差內(nèi)適用的最大水深為kh=5.1，表明文獻[2]公式也能滿足所有頻率的垂向積分誤差都低于1%.因此，針對D55工況，改進前后的計算結(jié)果不存在明顯差異.B55寬譜情況下，波數(shù)范圍在kh=1.568～7.825之間，結(jié)合圖1結(jié)果可知，在1%誤差內(nèi)，文獻[2]中水平速度計算公式不能涵蓋全部頻域，kh=7.825的誤差約為2%，本文水平速度在1%誤差內(nèi)適用的最大水深為kh=7.34，kh=7.825的誤差約為1.14%，可見改進公式的多數(shù)頻率范圍的誤差均能控制在1%左右，因此得到的速度剖面更為精確.在自由表面(波面)附近的結(jié)果比試驗結(jié)果大，一方面數(shù)值模擬中的衰減現(xiàn)象比物理模型試驗中的衰減現(xiàn)象弱；另一方面，物理模型試驗中的測量方法很難精確捕捉到自由表面處的速度，其他一些學者采用勢流理論模型進行求解也得到類似結(jié)果.在水底附近數(shù)值模擬的水平速度與試驗結(jié)果存在一定差異，主要是數(shù)值模擬采用Boussinesq水波方程，其在物理機制上與物理模型試驗中的衰減現(xiàn)象并不完全一致.

圖5 文獻[2]公式和改進公式的計算水平速度剖面與D55工況試驗結(jié)果的比較Fig.5 Comparisons of calculated horizontal velocity profiles from the present formula and the original formula with D55 experimental data

圖6 文獻[2]公式和改進公式的計算水平速度剖面與B55工況試驗結(jié)果的比較Fig.6 Comparisons of calculated horizontal velocity profiles from the present formula and the original formula with B55 experimental data

3 結(jié)論

提出一種改進方程速度精度的方法，即在雙層Boussinesq水波方程的速度公式基礎上，引入帶有常系數(shù)β的三階項，不改變原始Boussinesq水波方程，利用新的公式即可獲取更高精度的速度場.通過理論分析和數(shù)值驗證，主要得出以下結(jié)論：

(1) 雙層Boussinesq水波方程的速度場可以在一定程度上得到改進，常系數(shù)β=0.78時，在1%誤差內(nèi)，水平速度和垂向速度的最大適用水深得到了較大幅度改進，其中水平速度由原來的kh=5.1提高到kh=7.34，垂向速度由原來的kh=4.5提高到kh=7.83.

(2) 利用模型模擬強非線性聚焦波和流函數(shù)波浪(穩(wěn)態(tài)波)演化，發(fā)現(xiàn)改進的水平速度剖面與試驗結(jié)果的吻合度更高，從數(shù)值角度展示了理論改進方式的有效性.

此外，本文的改進方法可為其他Boussinesq水波方程改進提供重要參考，有關(guān)詳細的分析與討論有待更深入的研究.