基于非靜壓模型的波浪空間聚焦研究?

馬小舟， 馬建興， 馬玉祥， 董國海

(大連理工大學(xué)海岸和近海工程國家重點實驗室，遼寧 大連116024)

近年來極端波浪造成多起海洋事故，對海洋平臺以及海洋船舶的安全構(gòu)成極大威脅，而波浪聚焦是極端波浪的重要形成機(jī)理之一，因此對波浪聚焦的研究具有重要意義。波浪聚焦是指不同頻率、不同方向的波浪在傳播的過程中相互調(diào)制，使波能發(fā)生聚集，形成一個孤立的大波的過程。當(dāng)波浪經(jīng)過變水深區(qū)域時，如水底火山、海脊等，由于波浪折射作用波浪聚焦效果十分顯著。早期對該問題主要通過射線理論進(jìn)行研究[1]，但由于射線理論忽略了波浪繞射與衍射作用，在焦散線附近無法精確求解。對此，Berkhoff[2]在1972年推導(dǎo)出二維緩坡方程，綜合考慮了波浪折射和衍射現(xiàn)象，能顯著模擬波浪聚焦效果，計算聚焦區(qū)域波高?；谏鲜龌A(chǔ)，學(xué)者們針對不同地形的波浪聚焦問題開展了大量研究。

長波在海島附近的聚焦問題，Yu和Zhang[3]、Zhang和Zhu[4]、傅丹娟等[5]針對不同的海島簡化模型，研究發(fā)現(xiàn)在一定的條件下，島嶼附近會發(fā)生波浪“俘獲”現(xiàn)象，整個島嶼岸線都有較大波幅出現(xiàn)；Zhu和Harun[6]、Niu和Yu[7]研究了淺灘形狀、淺灘尺寸以及淺灘淹沒深度對淺灘附近波浪聚焦特性的影響；Longuet-Higgins[8]、Shaw和Neu[9]、熊夢婕等[10]對不同剖面海脊上俘獲波的解析解進(jìn)行了探討；Berkhoff等[11]、Vincent和Briggs[12]等通過物理實驗的方法對波浪通過淺灘后的波高變化進(jìn)行了研究。但現(xiàn)階段對于帶凸起的不規(guī)則斜坡地形的波浪聚焦問題少有研究。Whalin[13]在1972年用物理實驗的方法，在寬水槽中研究了凸起不規(guī)則斜坡地形上的波浪聚焦問題。但該實驗是在水槽中開展，由于水槽兩側(cè)存在波浪反射現(xiàn)象，因此并不是單純的波浪聚焦問題。本文在Whalin的實驗基礎(chǔ)上采用非靜壓模型對凸起不規(guī)則斜坡地形上的波浪聚焦特性進(jìn)行了研究。

1 模型簡介及驗證

本文選用由Zijlema等[14]建立的淺水非靜壓數(shù)值模型SWASH(Simulating Waves till SHore)對波浪聚焦問題進(jìn)行研究。該模型同時適用于波浪和水流條件，主要用于模擬沿岸地區(qū)的波浪和水流運動，是一個用于模擬二維或三維條件下的非靜壓、自由表面的旋流和波流輸運現(xiàn)象的數(shù)值工具。

1.1 控制方程及邊界條件

SWASH模型的控制方程采用的是常密度不可壓縮流體三維雷諾平均下的Navier-Stokes方程，在笛卡爾直角坐標(biāo)系(原點在靜水面上，z軸向上為正)中的張量表達(dá)式為：

(1)

(2)

(3)

(4)

Pt=Ph+Pd=g(ζ-z)+Pd。

(5)

方程在自由表面以及水底滿足如下方程：

(6)

計算過程忽略自由表面處風(fēng)應(yīng)力與張力，Pd|z=ζ=0，在底邊界處忽略兩個切向壓力。

入流邊界，由線性波浪理論得到波浪的法向速度分量，切向速度分量初始值為零。在出流邊界采用海綿層吸收波浪能以模擬開邊界，海綿層尺寸不計入計算域。在固邊界處，采用自由滑移的邊界條件。

1.2 數(shù)值求解

SWASH模型采用交錯網(wǎng)格有限差分法進(jìn)行空間離散，三維計算中速度變量定義在網(wǎng)格側(cè)面中心處，非靜壓 定義在網(wǎng)格上表面中心處，自由表面的壓力初始值為零，垂向分層一般選取1～3層。本文選用二階顯示有限差分法對水平動量方程做離散化處理，垂直動量方程選用Keller-box方法進(jìn)行離散化處理，在求解壓力速度耦合方程時采用分步法。垂向多層情況下的求解過程為由底邊界條件通過Keller-box離散方程以及連續(xù)方程由下向上求解垂向速度分量，然后將由水平動量方程、垂向速度以及自由表面動壓帶入連續(xù)方程通過Bi-CGSTAB求解下一時刻動壓，進(jìn)而求得下一時刻各個速度分量，最后根據(jù)由自由表面方程求解下一時刻的自由表面水位值 。更詳細(xì)的模型介紹可參見Stelling等[14]的相關(guān)文獻(xiàn)。

1.3 模型驗證

1.3.1 色散性驗證 SWASH模型中色散精度由垂向分層控制，垂向分層越多，模型性能越好，但因此造成的計算量也會大幅增長。本文中垂向分層均采用均勻分層。本節(jié)通過模擬規(guī)則波在二維平底地形上的傳播過程，驗證SWASH模型的色散精度。數(shù)值模擬過程水深取0.5 m，波高0.006 m，選取3種不同工況，使kd=π/4、π/3、π/2。沿波浪傳播方向網(wǎng)格步長取Δx=L0/50(L0為入射波波長)，垂向均勻分為3層。采用速度入口邊界，設(shè)定模擬過程中最大Courant數(shù)為0.5，水槽末端設(shè)置2倍波長的消浪區(qū)。距造波板10L0處三種不同工況下波浪模擬效果如圖1所示。圖中縱坐標(biāo)為無量綱化波面值，η為數(shù)值模擬波面，a0為初始波幅。由圖1可以看出，當(dāng)ka在0.003π～0.006π之間，kd取值在π/2～π/4之間時，SWASH模型垂向均勻分三層時可有效地模擬淺水以及中等水深中的波浪色散，效果理想。

圖1 數(shù)值模擬波面與理論波面對比(x=10L0)

1.3.2 非線性驗證 準(zhǔn)確模擬波浪的傳播除色散性以外另一個重要影響因素為非線性。本節(jié)采用SWASH模型模擬波浪在橢圓形淺灘地形上的傳播過程，通過比較淺灘附近的波高變化驗證SWASH模型的非線性模擬精度。實驗地形采用Vincent和Briggs[12]的實驗地形(見圖2)，橢圓形淺灘中心坐標(biāo)為x=6.10 m，y=13.72 m，地形及水深具體參數(shù)見Vincent和Briggs[12]。波浪由左側(cè)正向入射，入射波周期T=1.3 s，波高H=5.5 cm，計算域取25 m×30 m，由于波浪主要傳播方向為x方向，考慮到計算效率，計算網(wǎng)格沿x方向較密，y方向較稀疏，網(wǎng)格尺寸為0.05 m×0.1 m，計算區(qū)域末端設(shè)置5 m海綿層，用來吸收波浪能以模擬開邊界。垂向均勻分三層，模擬時長78 s。

圖2 實驗地形圖

為驗證SWASH模型模擬準(zhǔn)確性，選取波高變化較明顯的截面3、截面5和截面7上的波高分布與Vincent實驗數(shù)據(jù)進(jìn)行對比。通過上跨零點法統(tǒng)計數(shù)值模擬中截面處均方根波高值，與Vincent和Briggs[12]的實驗結(jié)果對比如圖3所示，圖中縱坐標(biāo)為均方根波高與入射波高的比值，實線為SWASH模擬結(jié)果，空心圓圈為Vincent實驗結(jié)果。由圖可以看出數(shù)值模擬結(jié)果與實驗結(jié)果吻合良好。

圖3 截面3、5、7上數(shù)模與Vincent實驗相對波高對比

2 數(shù)值模擬結(jié)果與討論

2.1 波浪聚焦驗證

在上述驗證基礎(chǔ)之上，為進(jìn)一步驗證SWASH軟件模擬波浪聚焦的準(zhǔn)確性，模擬Whalin[13]實驗過程。該實驗是在25.603 m×6.1 m的水槽(見圖4)中放置了一個三維聚焦地形，淺灘變水深區(qū)域等深線相互平行，在幾何上為直徑等于水槽寬度的半圓。波浪由深水處正向入射，水槽深水側(cè)水深為0.457 m，水槽末端淺水深為0.152 m。地形詳細(xì)表達(dá)式見文獻(xiàn)[13]。

圖4 實驗地形圖[13]

數(shù)值模擬中計算域取30 m×6.1 m，考慮波浪傳播方向以及計算效率，網(wǎng)格尺寸沿x方向較密，y方向較稀疏，為0.025 m×0.1 m，波浪從左側(cè)正向入射，入射波為規(guī)則波周期T=2 s，波高H=2.12 cm，在計算域末端設(shè)置10 m海綿層以吸收波浪能量。時間步長取0.005 s，垂向均勻分3層。

通過諧波分析的方法分析水槽中軸線上的波浪組成，將數(shù)值模擬結(jié)果與Whalin實驗結(jié)果進(jìn)行對比(見圖5)，其中橫坐標(biāo)為距造波板的距離，縱坐標(biāo)為無量綱化諧波波幅，實線為SWASH模擬結(jié)果，實心點為Whalin實驗結(jié)果。由圖5可以看出，數(shù)值結(jié)果與實驗結(jié)果吻合良好，隨著水深的減小，波浪在反射、折射以及淺水變形的作用下，基頻能量增大，聚焦區(qū)域由于波浪非線性二倍頻以及三倍頻能量迅速增大，聚焦區(qū)域后由于波浪折射作用波浪傳播方向發(fā)生改變，波浪迅速向水槽兩側(cè)傳播并與兩側(cè)壁的反射波相互作用波高減小。

圖5 數(shù)模與物模諧波分布圖對比

2.2 波浪聚焦特性的研究

波浪在斜坡地形上傳播時，當(dāng)?shù)匦蔚壬罹€與波浪傳播方向正交時，地形可簡化為二維地形，但當(dāng)斜坡起點不在同一條直線上波浪傳播方向與等深線斜交時，該地形則無法簡化成二維地形如Whalin實驗地形。但由于Whalin實驗中除波浪聚焦作用外還包含波浪反射的影響，為了研究單純的波浪聚焦特性，本節(jié)對Whalin實驗地形進(jìn)行調(diào)整如圖6所示。在原實驗地形的基礎(chǔ)上，將地形向兩側(cè)拓寬，該地形兩側(cè)斜坡起點橫坐標(biāo)x0=10.68 m，中軸線附近斜坡起點在水平面上為半徑等于3.05 m的半圓，借此消除波浪在水槽兩側(cè)的反射影響，波浪在該地形上的聚焦現(xiàn)象更接近實際地形上的波浪聚焦現(xiàn)象。

圖6 實驗地形圖

數(shù)值模擬中的入射波波浪要素以及地形水深條件均與Whalin實驗相同，波浪從左側(cè)正向入射，波浪入射側(cè)水深為0.457 m，淺灘頂水深為0.152 m，地形平行于x軸的截面均為1∶25的斜坡。計算域取35 m×18 m，網(wǎng)格尺寸為0.025 m×0.1 m，為吸收波浪能量在計算域兩側(cè)各設(shè)5 m海綿層，波浪入射方向設(shè)置10 m海綿層。數(shù)值模擬過程中的時間步長取0.005 s，垂向均勻分3層。

通過上跨零點法統(tǒng)計波浪在三維聚焦地形中軸線上的均方根波高，將其與入射波高取比值得到相對波高。相對波高沿水槽方向分布如圖7中黑色實線所示，藍(lán)色點劃線為Whalin實驗地形下中軸線上相對波高變化。由圖可以看出，水槽兩側(cè)的波浪反射現(xiàn)象不僅導(dǎo)致計算域內(nèi)波面分布現(xiàn)象不同，波高變化量也不同，反射波減小了聚焦區(qū)域內(nèi)的相對波高。同時由于波浪聚焦過程中包含淺水變形的影響，因此將相同入射波浪在同等水深變化的二維地形上的波高變化用紅色虛線表示。由圖可以看出波浪在二維地形上傳播時，由于淺水變形相對波高在淺灘頂出現(xiàn)最大值1.34H0(H0為入射波高)，波浪在三維聚焦地形上傳播時由于水深變化以及波浪折射作用波能聚集，最大波高出現(xiàn)在二維淺水變形中的最大波高之后，其值為2.81H0，而由于反射波作用的Whalin實驗地形上相對波高變化介于兩者之間。由結(jié)果可以看出，波浪在有凸起的斜坡淺灘上傳播時，波高變化中只有一小部分由于淺水變形導(dǎo)致，而波高變化的最主要原因為波能聚集。也就是說有凸起的斜坡淺灘上波浪聚焦是波高變化的最主要原因。

圖7 不同地形下相對波高對比

為進(jìn)一步研究波浪聚焦特性，對二維地形以及三維聚焦地形中軸線上的波浪進(jìn)行諧波分析如圖8所示，圖中黑色實線為相對波高的沿程分布，彩色虛線為一倍頻、二倍頻、三倍頻以及四倍頻的沿程分布，黑色豎線表示最大相對波高所在位置。通過諧波分析的結(jié)果可以看出，波浪在二維地形上傳播時，隨著水深變化，基頻、二倍頻、三倍頻以及四倍頻能量依次開始增大，隨著波浪傳播非線性使諧波之間存在相互能量交換，基頻能量減小，倍頻能能量增加。二倍頻、三倍頻以及四倍頻同時達(dá)到最大值，同時基頻達(dá)到極小值，此后基頻與倍頻均在同一位置達(dá)到極值。最大相對波高位于基頻最大值以及倍頻最大值之間。波浪在三維聚焦地形上傳播時，二倍頻、三倍頻以及四倍頻同樣在同一位置達(dá)到最大值，但基頻最大值位于倍頻最大值之后。與波浪在二維地形上傳播時相同，非線性同樣使諧波之間存在相互能量交換，但基頻極值與倍頻極值出現(xiàn)在不同位置。最大相對波高位于倍頻最大值與基頻最大值之間。通過比較可以看出波浪聚焦主要影響波浪的基頻能量及其分布，倍頻能量相比二維地形上大小變化明顯，但極值位置變化不大。

圖8 諧波分布圖

2.3 波浪聚焦的影響因素

通過上節(jié)討論可知，波浪在有凸起的斜坡淺灘上波浪聚焦現(xiàn)象十分顯著。本節(jié)通過SWASH模型模擬同波陡不同波長的波浪在該三維聚焦地形上的傳播，進(jìn)一步研究初始無量綱參數(shù)kR對波浪聚焦位置的影響。本文中選取較小的入射波陡值保證計算域內(nèi)出現(xiàn)的最大波浪未達(dá)到破碎指標(biāo)，計算域內(nèi)無波浪破碎。本節(jié)中共設(shè)置了7種同波陡不同波長的波況，波浪參數(shù)如表1所示。實驗參數(shù)設(shè)置與2.2節(jié)中三維聚焦地形相同。

表1 波浪參數(shù)

數(shù)值模擬過程中，通過對波浪聚焦模擬過程的動畫以及剖面波高進(jìn)行觀測，確保波浪聚焦達(dá)到的最大波高未發(fā)生破碎。對中軸線上最大波高位置處的波形進(jìn)行監(jiān)測，如圖9所示。圖中分別為工況A01、A03、 A06下中軸線上最大波高位置處的波浪形態(tài)，其中橫坐標(biāo)為無量綱時間，縱坐標(biāo)為無量綱化波面，其中a0為入射波波幅，T為入射波浪周期。由圖可以看出隨著初始kR的減小，最大波高處波浪不對稱性增大，非線性增強，聚焦區(qū)域距造波板越來越近。當(dāng)初始kR=1.56π時，波浪在非線性作用下波浪形態(tài)由單峰變?yōu)殡p峰。

圖9 工況A01、A03、A06最大波高處波浪形態(tài)

通過上跨零點法統(tǒng)計不同波況下計算域內(nèi)的均方根波高，與入射波高取比值得到相對波高。相對波高等值線在聚焦地形附近的分布情況如圖10所示。圖中分別為波況A03和A06的相對波高等值線分布情況。對不同工況下最大相對波高距造波板的距離進(jìn)行統(tǒng)計，最大相對波高位置與初始kR的關(guān)系如圖11所示，圖中縱坐標(biāo)表示最大相對波高距造波板的距離，橫坐標(biāo)為初始kR值。由圖10、11可以看出初始無量綱參數(shù)kR越大，聚焦區(qū)域越小，同時最大波高的位置距造波板越遠(yuǎn)，這是由于波浪在等深線為圓形的地形上傳播時，入射角相同時初始kR越小，折射作用下波向線偏轉(zhuǎn)速度越快[15]，因此聚焦區(qū)域距聚焦地形越近。

圖10 聚焦區(qū)域相對波高等值線分布圖

圖11 最大相對波高位置與初始波數(shù)的關(guān)系

由前面的分析可知，初始kR越大聚焦區(qū)域越小，如果僅考慮聚焦區(qū)域的大小對相對波高的影響，則聚焦區(qū)域越小平均相對波高越大。但波浪能量在聚焦區(qū)域并非均勻分布因此不能僅僅通過聚焦區(qū)域的大小確定最大相對波高的大小，還需要考慮波浪能量在聚焦區(qū)域的分布以及不同波浪間的相互作用。通過圖9可以看出線性波浪隨著非線性增加，波浪形態(tài)逐漸偏離正弦分布，波峰變尖波谷變平坦，但當(dāng)非線性進(jìn)一步增大，倍頻能量增加時波浪將由單峰變?yōu)殡p峰甚至多峰，因此在最大波浪未達(dá)到破碎條件的情況下，最大相對波高的大小除需要考慮聚焦區(qū)域大小外，還要考慮波浪非線性的影響。圖12給出了不同工況下最大相對波高的值與初始kR的關(guān)系。由圖可以看出當(dāng)初始kR在1.4π～4.05π之間時，隨著初始kR的減小，最大相對波高先增大后減小，在kR=2.45π時，最大相對波高達(dá)到極大值，可以達(dá)到2.48H0(H0為入射波波高)。

圖12 最大相對波高與初始波數(shù)的關(guān)系

3 結(jié)論

(1) 非靜壓模型SWASH軟件通過將垂向壓力分解為動壓和靜壓以及區(qū)別于其他非靜壓軟件的動壓位置定義，采用三層垂向分層即可保證足夠的色散精度以及非線性精度，準(zhǔn)確高效的模擬波浪在淺水以及有限水深中的傳播。

(2) 通過SWASH模型模擬波浪在三維聚焦地形上的傳播，將數(shù)值模擬結(jié)果與物理實驗結(jié)果進(jìn)行對比，證明了SWASH模型在模擬地形導(dǎo)致的波浪折射、反射、衍射以及淺水變形等多重作用下的準(zhǔn)確性。

(3) 采用SWASH模型模擬相同波浪在二維斜坡地形以及三維凸起斜坡淺灘地形上的傳播，分析了相對波高的沿程分布以及諧波組成，對比發(fā)現(xiàn)，波浪非線性使諧波之間存在相互能量交換，波浪在二維地形上傳播時，波浪傳播過程中波能守恒，波浪在三維地形上傳播時，波浪聚焦波浪能量增加，非線性增強，波浪聚焦引起的波高變化遠(yuǎn)大于二維淺水變形引起的波高變化。

(4) 通過SWASH模型模擬同波陡不同波長的波浪在相同三維地形上的傳播，研究了最大波浪無破碎的情況下，聚焦區(qū)域的位置與初始無綱參數(shù)kR的關(guān)系，以及最大相對波高的影響因素。結(jié)果表明，在最大浪未發(fā)生破碎的情況下，等深線為圓形的聚焦地形下初始無綱參數(shù)kR對波浪聚焦有顯著影響，初始kR值越大，波向線偏轉(zhuǎn)速度小，最大波高位置距斜坡起點距離越大。波浪聚焦產(chǎn)生的最大波高值不僅與聚焦區(qū)域大小有關(guān)，還與波浪非線性的強弱有關(guān)，當(dāng)kR在1.4π～4.05π之間時，隨著無綱參數(shù)kR的減小，最大相對波高先增大后減小。kR=2.45π時，最大相對波高出現(xiàn)極大值，可以達(dá)到2.48H0。

(5) 本文通過SWASH模型對規(guī)則波在半圓形凸起斜坡地形上的聚焦特性進(jìn)行了研究，但未考慮最大波浪發(fā)生破碎對波浪聚焦特性的影響。未來還需對波浪破碎對波浪聚焦的影響進(jìn)行系統(tǒng)研究。并且不規(guī)則波在凸起斜坡地形的聚焦特性，斜坡凸起的形狀是否會對聚焦特性有所影響？波浪聚焦過程的能量時空變化是怎樣，這些問題需要進(jìn)一步深入分析與討論。