崔　成，肖　輝，左書華，張義豐

(1.交通運輸部 天津水運工程科學研究所，天津 300456；2. 港口水工建筑技術國家工程實驗室，天津 300456)

?

數(shù)值研究非平底地形對畸形波能量結構的影響*

崔成1,2，肖輝1,2，左書華1,2，張義豐1,2

(1.交通運輸部 天津水運工程科學研究所，天津 300456；2. 港口水工建筑技術國家工程實驗室，天津 300456)

采用流體體積(VOF)方法捕捉自由表面，結合有限差分方法求解N-S(Navier-Stokes)方程、k-ε模型封閉，建立數(shù)值模型，并使用該模型模擬非平底地形條件下畸形波的生成。采用小波分析方法計算模擬結果的時頻能量譜，基于計算結果分析非平底地形對畸形波能量集中程度和高頻能量的影響。主要結論為：坡度小于1∶10的斜坡地形和無量綱水深變化小于0.333的曲線地形對畸形波時頻能量譜的影響不顯著；坡度大于1∶10的斜坡地形和無量綱水深變化大于0.333的曲線地形會顯著影響畸形波的時頻能量譜，隨著坡度和無量綱水深變化的增加，時頻能量譜中畸形波發(fā)生時刻附近，能量向高頻方向移動，使得能量在高頻端的分布范圍增大，時頻譜密度峰值減??；斜坡和曲線地形的特征變化對于畸形波能量集中度參數(shù)的影響不顯著。

畸形波；能量結構；非平底地形

畸形波是一種包含在隨機波列當中的單個異常大波，具有明顯的非線性特征，持續(xù)時間短，能量集中，破壞力驚人[1]。隨著人類對海洋的開發(fā)和利用不斷發(fā)展，關于畸形波的研究已經引起越來越多國內、外學者的關注。目前，畸形波的生成機理還沒有定論，根據(jù)畸形波發(fā)生的環(huán)境條件以及內部、外部特征，學者們提出了一些假說，如：組成波波能聚焦作用、波流相互作用、大氣與海水的相互作用、地形作用、波浪調制不穩(wěn)定機制、波浪間的非線性相互作用等[2]，這些假說在逐漸被證明的過程中。

上述可能引發(fā)畸形波的因素各自并非獨立存在，而且之間也存在復雜的聯(lián)系。因此一次畸形波事件可能不僅是一種生成機制引發(fā)，可能為上述兩種或兩種以上生成機制。對于近岸較淺水域，地形變化條件普遍存在，如果滿足有限水深條件(水深與特征波長之比小于0.5)，則地形對波浪傳播、演變的影響則不能忽視。當波浪經過復雜地形從較深水域傳播至較淺水域的過程中，復雜地形對波浪的反射，折射和變淺作用會改變波浪的傳播形態(tài)和能量傳遞過程。Biausser等[3]和Guyenne等[4]數(shù)值模擬了孤立波在水深不斷變淺的非平底地形上的傳播演化過程，結果顯示未發(fā)生破碎之前，隨著水深的不斷減小，能量不斷的積累，波高和波陡不斷的增加。Sergeeva等[5]使用KDV方程模擬了隨機波列在不規(guī)則地形上的傳播和變形，研究結果發(fā)現(xiàn)波高隨地形的起伏而變化，且波高的分布偏離Rayleigh分布，大波的發(fā)生概率隨地形突變個數(shù)的增多而增大；相對于平底地形，當波列在不規(guī)則地形上傳播演化時，畸形波的發(fā)生概率明顯提高。裴玉國[6]的物理試驗結果也得到了相似的結論。由此可見地形作用可能是引發(fā)畸形波的因素之一。在有限水深條件下，研究水底地形對畸形波的影響有助于解釋畸形波的生成機理?；尾ǖ膬炔磕芰拷Y構決定其外部特征，能更加直觀地描述畸形波的形成機制。小波分析方法可以有效地捕捉畸形波發(fā)生時的能量分布狀態(tài)，已有學者使用小波分析方法開展了有關畸形波時頻能量譜的研究。主要研究結果表明，畸形波的小波譜具有十分明顯的特征，能量集中且包含大量的高頻成分，非線性相互作用明顯，部分實測畸形波的小波譜還表現(xiàn)出雙峰特征[7-11]。由此可見，對畸形波能量結構的研究可以重點關注高頻能量成分和能量集中程度兩個方面。

本文采用流體體積(VOF)方法捕捉自由表面，并結合有限差分方法求解N-S(Navier-Stokes)方程建立數(shù)值模型；運用雙波列疊加模型(組成波聚焦機制)模擬非平底地形條件下(地形作用)畸形波的生成；采用小波分析方法計算模擬畸形波的時-頻能量譜；根據(jù)計算結果進一步分析地形條件變化對其內部能量特征的影響。

1數(shù)學模型建立和驗證

1.1控制方程

考慮到畸形波的波陡很大，在其形成過程中能量不斷匯聚，流體湍動效應顯著，因此考慮適當?shù)耐牧髂Ｐ汀Ψ欠€(wěn)態(tài)的Navier-Stokes方程進行時間平均運算，得到雷諾時均方程，并利用目前應用范圍廣、檢驗程度高的k-ε模型建立紊流脈動值附加項與其他時均值之間的聯(lián)系，來封閉雷諾時均方程，與連續(xù)方程一起作為控制方程，見式(1)～(5)。

連續(xù)方程：

(1)

雷諾時均方程：

(2)

(3)

雙方程k-ε模型：

(4)

(5)

式中，u和v分別表示x和y方向的速度分量；θ表示部分單元參數(shù)，不依賴于時間，依據(jù)網(wǎng)格中流體的比例，取值為0～1；p為流體壓力；ρ為流體的密度；ν為流體的運動黏滯系數(shù)；νt=Cu(k2/ε)，為湍動黏滯系數(shù)；k表示湍動能；ε表示湍動能耗散率。k-ε方程中的經驗參數(shù)分別選為：Cu=0.09，σk=1.0，σε=1.3，σε1=1.43，Cε2=1.92[12]。

1.2邊界條件

1.2.1造波邊界條件

造波邊界條件采用雙波列疊加模型計算，見式(6)。該模型中，目標譜能量被分成2部分：一部分能量分配給瞬態(tài)波列，另一部分分配給常規(guī)隨機波列[13]。

(6)

式中，η(x,t) 表示距離造波板x處的波面時間過程；η1(x,t) 和η2(x,t) 分別表示常規(guī)隨機波列和瞬態(tài)波列的波面時間過程；ki，ωi(ωi在第i個頻域區(qū)間上隨機選取，為避免波浪以周期2π/Δω重復出現(xiàn))和εi(0～2π內均勻分布)分別表示第i個組成波的波數(shù)，角頻率和隨機初相位；xc和tc分別表示瞬態(tài)波的聚焦時間和位置；a1i和a2i分別表示兩波列第i個組成波的振幅，由式(7)確定：

(7)

式中，p1=80%，p2=20%，該分配比例能生成滿足定義的畸形波[13]；S(f)表示譜密度函數(shù)，本文選用修正的P-M譜[14]作為目標譜，表達式見式(8)：

(8)

式中，T表示平均周期;Hs表示有效波高。

根據(jù)造波原理，造波邊界的水平速度U0(t)可表示為

(9)

1.2.2吸收邊界條件

為了限制計算域，在水槽的右端設置一個吸收邊界層來吸收出流波浪。在吸收邊界區(qū)域中，速度和自由表面等變量分別衰減為：v(x,y)/μ(x)，η/μ(x)，k/μ(x)和ε/μ(x)。μ(x)為衰減系數(shù)：

(10)

式中，λ為吸收區(qū)域的長度；λ′為吸收區(qū)域內的點到水槽右端邊界的距離；α為阻尼系數(shù)。

1.3數(shù)值模型驗證

為了驗證所建立數(shù)值模型模擬畸形波的有效性，使用該數(shù)值模型模擬一組畸形波。計算區(qū)域如圖1所示，在x和y方向分別劃分1 600和30個單元格，計算域右端設置30 m吸收區(qū)域。在x和y方向單元格間距分別為5和3 cm，計算時間220 s，時間間隔0.01 s，運動黏性系數(shù)1.002×10-6m2/s。驗證采樣點距離造波邊界2 100 cm。數(shù)值模擬中畸形波模擬方法及驗證點位置均與物理試驗相匹配。物理試驗在大連理工大學海岸及近海工程國家重點實驗室的渾水水槽中進行，該水槽長56 m，寬3 m，最大工作水深0.7 m。圖2給出了數(shù)值模擬畸形波波面時間過程與物理實驗結果的對比。從圖2的對比可以看出，本文所建立的數(shù)值模型可以較準確地模擬具有較強非線性的畸形波。圖2b為局部放大圖。

圖1　數(shù)值模型計算域的布置圖(cm)Fig.1　Setup of the numerical modeling domain

圖2　畸形波波面時間過程數(shù)值模擬與物理實驗結果的對比Fig.2　Comparison of the simulated and experimental freak waves

為了驗證該數(shù)值模型模擬非線性波浪在非平底地形上傳播的有效性。使用該數(shù)值模型模擬孤立波在1∶20斜坡上傳播的分裂現(xiàn)象，水深從d0=7.62 cm 逐漸變?yōu)閐1=3.81 cm。計算域與Madsen和Mei[15]所做的物理實驗相匹配，見圖3。孤立波波要素為：振幅0.91 cm，傳播速度91.96 cm/s。邊界條件，網(wǎng)格尺寸和時間步長都與Choi和Wu[16]給出的數(shù)值模擬算例一致。圖4給出了模擬結果對比，為了方便對比統(tǒng)一將a位置波峰對應時刻設置為0 s?？傮w看來，三者基本一致，在水深逐漸變淺的斜坡上，波面發(fā)生了分裂現(xiàn)象，在孤立波的右側出現(xiàn)了一個較小的波浪。與Choi和Wu[16]數(shù)值結果的差異是由于本文數(shù)值模型考慮黏性耗散的影響，而前者未考慮黏性耗散。比較而言，本文的模擬結果無論是相位還是峰值，與Choi和Wu的數(shù)值結果相比都更為接近物理模擬結果。上述結果表明，本文所建立的數(shù)值模型模擬較強非線性波浪在地形上傳播和演化過程具有較好的精度和穩(wěn)定性。

圖3　孤立波在斜坡地形上傳播的布置圖(cm)Fig.3　Setup of a solitary wave propagation over an non-flat bottom

圖4　孤立波傳播、演化過程的數(shù)值結果驗證示例Fig.4　Comparison of the modeled and the experimental solitary wave propagating over a non-flat bottom

2　數(shù)值計算組別

本文考慮兩種非平底地形，分別為斜坡地形和拋物曲線地形。對于斜坡地形，模擬了8組不同坡度的斜坡，坡度在1∶25 ～ 1∶5.6范圍內變化。斜坡地形邊界和計算域參見圖5。將數(shù)值水槽分為3部分，設造波邊界為原點，X=0～1 250 cm為平底地形，水深為d0=150 cm；X=1 250～2 000 cm為斜坡地形，水深從d0變化到d1，對應8組坡度坡頂?shù)乃頳1在120～15 cm變化；X=2 000～6 500 cm為平底地形，水深為d1。對于拋物曲線地形，模擬了6組拋物曲線地形，地形最大高度在30～130 cm范圍內變化。拋物曲線地形邊界和計算域參見圖6。X=0～1 250 cm為平底地形，水深為d0=150 cm；X=1 250～2 750 cm為拋物曲線地形，水深從d0變化到d1，對應6組曲線地形，坡頂?shù)乃頳1在120～20 cm變化，地形曲線由拋物線方程確定；X=2 750～6 500 cm為平底地形，水深為d0。

圖5　斜坡地形的計算域示意圖Fig.5　Sketch of the numerical modeling domain with a slope

圖6　曲線地形的計算域示意圖(cm)Fig.6　Sketch of the modeling domain with bump bottom

表1匯總給出了數(shù)值模擬算例組別。表中Ls為有效波長(有效周期對應的波長)；Δd/d0為無量綱水深變化，其中Δd=d0-d1。表2給出模擬畸形波的特征參數(shù)。

所有模擬組別中，入射波條件相同，均為有效波高Hs=5 cm，有效周期Ts=1.4 s。無論平底地形(Case 1)還是非平底地形(Case 2～Case 15)均得到了滿足定義的畸形波，即波高H0大于等于波列有效波高Hs的2倍，α1=H0/Hs≥2，波高H0大于等于與該波浪相鄰波浪波高的2倍，α2=H0/H-1≥2和α3=H0/H+1≥2，波峰η0大于等于波高的0.65倍，α4=η0/H0≥0.65。

表1　數(shù)值模擬組別

表2　畸形波特征參數(shù)

3　小波譜計算和分析

3.1小波譜計算方法

將包含畸形波的連續(xù)波浪序列做連續(xù)小波變化：

(11)

(12)

定義無量綱時頻能量譜密度S*(f*,T*)：

(13)

式中，f*為無量綱頻率，f*=f/fp；f表示頻率；fp表示譜峰頻率；T*為無量綱時間；T*=(t-tc)/Tp；tc表示畸形波生成、演化過程中發(fā)生異常大波波峰對應時刻，Tp表示譜峰周期。

3.2小波譜分析

圖7　平底和斜坡地形條件下畸形波的無量綱時頻能量譜Fig.7　Non-dimensional time-frequency spectrums of freak waves for cases 1～9

圖7給出了平底地形(Case 1)及8組斜坡地形(Case 2～Case 9)條件下畸形波的時頻能量譜。從圖7中可以看出：當坡度s從0增長到1∶10時，畸形波的時頻能量譜基本沒有變化；當坡度s從1∶10增長到1∶5.6時，在畸形波發(fā)生時刻附近，能量有向高頻端移動的趨勢；無量綱時頻譜密度峰值由30～35降到20～30；能量在高頻端的分布范圍變大，在坡度為1∶5.6的斜坡地形條件下，高頻端能量分布范圍可達3倍的譜峰頻率(平底地形條件下，高頻端能量分布范圍為2.5倍的譜峰頻率)。

圖8給出了平底地形(Case 1)及6組拋物曲線地形(Case 10～Case 15)條件下畸形波的時頻能量譜。由圖8可知：當Δd/d0從0增長到0.333時，畸形波的時頻能量譜基本沒有變化；當Δd/d0從0.333增長到0.867時，在畸形波發(fā)生時刻附近，能量有向高頻端移動的趨勢；無量綱時頻譜密度峰值由30～35降到20～30；能量在高頻端的分布范圍變大，在Δd/d0=0.867的曲線地形條件下，高頻端能量的分布范圍可達4倍的譜峰頻率。

圖8　平底和曲線地形條件下畸形波的無量綱時頻能量譜Fig.8　Non-dimensional time-frequency spectrums of freak waves for cases 1, 10～15

由此可見，坡度s小于1∶10的斜坡地形和無量綱水深變化Δd/d0小于0.333的曲線地形，對畸形波時頻能量譜的影響不顯著；而對于坡度大于1∶10的斜坡地形和無量綱水深變化Δd/d0大于0.333的曲線地形，會顯著影響到畸形波的時頻能量譜。隨著s和Δd/d0的增加，在畸形波時頻能量譜中，畸形波發(fā)生時刻附近，能量向高頻端移動，使得能量在高頻端的分布范圍增大，時頻譜密度峰值減小。

為了進一步量化分析地形變化對畸形波時頻能量譜的影響，圖9給出了畸形波能量集中程度參數(shù)αE隨地形特征變化的情況。

其中畸形波能量集中程度參數(shù)αE定義為對于某特定時刻，各頻率成分的總能量Ec：

(15)

圖9　畸形波αE隨坡度的變化Fig.9　Relationship between the slope s and αE of freak waves for cases 1～15

(16)

則能量集中度參數(shù)αE：

(17)

從圖9中可以看出，隨著Δd/d0的增加，畸形波能量集中度參數(shù)αE變化很小，兩種地形條件下能量集中度參數(shù)αE的變化情況一致。

除能量集中度外，另一畸形波能量特征為包含較多高頻成分，為了定量分析地形變化對畸形波高頻能量的影響，計算15組算例畸形波周期時段內，2倍譜峰頻率以上頻率成分能量和Eh及各頻率成分總能量Et。圖10給出了2種地形條件下，畸形波Eh/Et隨地形特征的變化,其中s表示斜坡坡度。

圖10　畸形波Eh/Et隨地形特征的變化Fig.10　Relationship between the slope s and Eh/Et of the freak waves for cases 1～15

從圖10中可以看出：斜坡地形條件下s≤1∶6.5時和曲線地形條件下Δd/d0≤0.667時，地形變化對Eh/Et的影響不顯著(相對于平底地形算例的變化幅度≤10%)；對于斜坡地形條件，當坡度s≥1∶6.0時，變化幅度有一定程度的增長，在20%～40%之間，Eh/Et隨著坡度的加大而增加；對于曲線地形條件，當Δd/d0≥0.833時，變化幅度增長明顯，在90%～310%之間，Eh/Et隨著坡度的加大而增加。曲線地形條件下無量綱水深變化Δd/d0對Eh/Et的影響遠遠大于斜坡地形條件下坡度s對Eh/Et的影響。

4　結　論

基于數(shù)值模擬結果，采用小波分析方法計算畸形波的時-頻能量譜，通過分析斜坡地形和曲線地形特征變化對時-頻能量譜的影響主要得到如下結論：

1)坡度小于1∶10的斜坡地形和無量綱水深變化小于0.333的曲線地形，對畸形波時頻能量譜的影響不顯著；對于坡度大于1∶10的斜坡地形和無量綱水深變化大于0.333的曲線地形，對畸形波的時頻能量譜有一定影響。隨著坡度和無量綱水深變化的增加，畸形波發(fā)生時刻附近，能量向高頻端移動，使得能量在高頻端的分布范圍變大，時頻譜密度峰值減小。

2)斜坡和曲線地形的特征變化對于畸形波能量集中度參數(shù)αE的影響不顯著。

3)坡度小于1∶6.5和無量綱水深變化小于0.667算例，地形變化對2倍譜峰頻率以上高頻成分的影響不顯著；當坡度大于1∶6.0和無量綱水深變化大于0.833時，隨著坡度和無量綱水深變化的增加，2倍譜峰頻率以上高頻成分明顯增加。

[1]CUI C. Study on generation、evolution and internal structures of freak waves[D].Dalian:Dalian University of Technology, 2013.崔成.畸形波生成、演化及內部結構研究[D].大連:大連理工大學,2013.

[2]KHARIF C, PELINOVSKY E. Physical mechanisms of the rogue wave phenomenon[J].European Journal of Mechanics B/Fluids, 2003,22:603-634.

[3]BIAUSSER B, GRILLI S T, FRAUNIE P. Numerical analysis of the internal kinematics and dynamics of three-dimensional breaking waves on slopes[C]∥Proceeding 13th Offshore and Polar Engineering Conference, USA, 2003:340-346.

[4]GUYENNE P, GRILLI S T. Computations of 3D overturning waves in shallow water[C]∥Proceeding 13th Offshore and Polar Engineering Conference, USA, 2003:347-352.

[5]SERGEEVA A, PELINOVSKY E, TALIPOVA T. Nonlinear random wave field in shallow water: variable Korteweg-de Vries framework[J].Natural Hazards and Earth System, 2011,11(2):323-330.

[6]PEI Y G. The generation of freak waves and its Behaviors[D].Dalian: Dalian University of Technology, 2007.裴玉國.畸形波的生成及基本特性研究[D].大連:大連理工大學,2007.

[7]LIU P C, MORI N. Wavelet spectrum of freak waves in the ocean[C]∥Proceedings, 27th International conference on Coastal Engineering, ASCE, 2000:1092-1098.

[8]CHIEN H, KAO C C. On the characteristics of observed coastal freak waves[J].Coastal Engineering Journal, 2002,44(4):301-319.

[9]CHRISTOU M., EWANS K., BUCHNER B, et al. Spectral characteristics of an extreme crest measured in a laboratory basin[J/OL]. http:∥www.ifremer.fr/web-com/stw2008/rw/full/Christou.pdf.

[10]ZHAO X Z. Experimental and numerical study of freak waves[D].Dalian: Dalian University of Technology, 2008. 趙西增.畸形波的實驗研究和數(shù)值模擬[D].大連:大連理工大學,2008.

[11]CUI C, ZHANG N C. Research on the time-frequency energy structure of generation and evolution of freak wave[J].The Ocean Engineering, 2011,29(3):59-66.崔成,張寧川.畸形波生成、演化過程時頻能量結構研究[J].海洋工程,2011,29(3):59-66.

[12]RODI W. Turbulence models and their application in hydraulics[C]∥IAHR Monograph，1984(8)：58.

[13]KRIEBEL D L. Efficient simulation of extreme waves in a Random Sea[J/OL].http:∥www.ifremer.fr/metocean/conferences/stw_abstracts/kriebel.pdf.

[14]YU Y X, LIU S S. Random wave and its applications to engineering[M].Dalian: Dalian University of Technology Press, 2011.俞聿修,柳淑學.隨機波浪及其工程應用[M].大連:大連理工大學出版社,2011.

[15]MADSEN O S, MEI C C. The transformation of a solitary wave over uneven bottom[J].Journal of Fluid Mechanics, 1969,39:781-791.

[16]CHOI D Y, WU C H. A new efficient 3D non-hydrostatic free-surface flow model for simulating water wave motions[J].Ocean Engineering, 2006,33:587-609.

[17]TORRENCE C, COMPO G P. A practical guide to wavelet analysis[J].Bulletin of the American Meteorological Society, 1998,79(1):61-78.

Numerical study on the effects of non-flat bottom topography on the energy structure of freak waves

CUI Cheng1,2，XIAO Hui1,2，ZUO Shu-hua1,2，ZHANG Yi-feng1,2

(1.TianjinResearchInstituteforWaterTransportEngineering，M.O.T.,Tianjin 300456,China;2.NationalEngineeringLaboratoryforPortHydraulicConstructionTechnology，Tianjin 300456，China)

A numerical model that solves the Reynolds time-averaged N-S equations by the finite-difference method was built with the k-ε model for the turbulence closure and the volume of fluid (VOF) method for tracking the free surface. The model was implemented to simulate the freak waves over non-flat bottom topography. With the time-frequency energy spectrum of the freak waves that was calculated by the wavelet analysis method the effects of non-flat bottom topography on the energy focusing and high frequency energy were analyzed. It is concluded that the non-flat bottom topography has a insignificant effect on the time-frequency energy spectrums for the slope that is less than 1∶10 and the topography whose normalized change in water depth is smaller than 0.333, but the influence of the non-flat bottom topography is considerable for the slope that is more than 1∶10 and the topography whose normalized change in water depth is greater than 0.333. With the increasing slope and normalized change in water depth, the energy shifts to high frequency endwhich makes the frequency domain of the energy distribution broaden, and the peak value of the time-frequency spectrum density decreases. On the other hand, the non-flat bottom topography has a indifferent effect on the parameters of energy focusing.

freak wave; energy structure; non-flat bottom topography

August 18, 2015

2015-08-18

國家自然科學基金項目——量化分析三維畸形波內部結構(51509120)

崔成(1984-)，男，遼寧海城人，博士，主要從事波浪、潮流模擬及其與結構物相互作用方面研究.E-mail：chengcui1984@163.com

(陳靖編輯)

P731

A

1671-6647(2016)03-0377-11

10.3969/j.issn.1671-6647.2016.03.007