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

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

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

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

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

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

1 速度公式的修正

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

從靜水位到連接面的速度場(chǎng)為

(1)

(2)

從連接面到水底的速度場(chǎng)為

(3)

(4)

從自由表面到靜水位的速度場(chǎng)為

(5)

(6)

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

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

(7)

(8)

(9)

(10)

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

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

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

(11)

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

(12)

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

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

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

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

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

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

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

2.1 數(shù)值模型

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

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