鄒 麗 , 孫 澤 , 宗 智 , 王 振 , 蔡志文 , 劉小龍
(1.大連理工大學(xué) 船舶工程學(xué)院,遼寧 大連116024;2.中國船舶科學(xué)研究中心,江蘇 無錫214082;3.工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室,遼寧 大連 116024;4.大連理工大學(xué) 數(shù)學(xué)科學(xué)學(xué)院,遼寧 大連116024)
波浪在復(fù)雜海洋地形中傳播時(shí),需要經(jīng)過多變的海底和復(fù)雜島礁地形;而島礁附近一般有幾公里的珊瑚礁淺灘,這導(dǎo)致流和波浪在島礁附近的繞流和傳播極其復(fù)雜,存在過頂繞流和三維繞流現(xiàn)象,伴隨有爬坡、折射,繞射、波浪破碎等非線性復(fù)雜問題。對(duì)海洋島礁近岸水波的數(shù)學(xué)建模應(yīng)該能夠準(zhǔn)確地描述水波運(yùn)動(dòng)的主要特征(波形的傳播,色散,波的折射、反射和繞射,波能的傳播,淺水效應(yīng),波的衰減,擴(kuò)散,破碎、爬坡等)。
基于能量平衡的波浪計(jì)算模型,或稱為波能譜或波作用密度譜模型。該模型基于波浪傳播過程中波譜作用量平衡原理,建立求解方程。它主要用于深海波浪的計(jì)算,亦適用于近海大范圍的波浪計(jì)算。本文針對(duì)島礁周邊遠(yuǎn)場(chǎng)較深水深的波浪傳播,采用動(dòng)譜平衡方程模型中的WAVEWATCH-III模型處理。WAVEWATCH-III是在NOAA/NCEP發(fā)展的第三代波浪模型,也是基于NASA戈達(dá)德宇宙飛行中心研發(fā)出的WAVEWATCH-II,代爾夫特理工大學(xué)研發(fā)的WAVEWATCH-I模型發(fā)展而來的。WAVEWATCH-III不僅僅是一個(gè)波浪模型,它更像是一個(gè)處理波浪相關(guān)問題的平臺(tái),科研工作者允許自行開發(fā)解決波浪問題的物理和數(shù)值方法。許多學(xué)者利用風(fēng)場(chǎng)資料驅(qū)動(dòng),在WAVEWATCH-III模型下對(duì)各大海域進(jìn)行的數(shù)值模擬預(yù)報(bào),取得了符合實(shí)測(cè)資料的結(jié)果。Feng和Vandemark等(2004)[1]將該模型應(yīng)用于修正衛(wèi)星測(cè)量的海平面高度,分析對(duì)比了純NCEP風(fēng)場(chǎng)和NCEP/QSCAT混合風(fēng)場(chǎng)驅(qū)動(dòng)下模型計(jì)算所得的波浪場(chǎng),包括了全球范圍和局部海洋范圍內(nèi)的情況。Browne等(2006)[2]用WAVEWATCH-III結(jié)合人工神經(jīng)網(wǎng)絡(luò)的經(jīng)驗(yàn)方法估算了近岸的波浪,認(rèn)為這種方法比直接應(yīng)用近岸模型SWAN的效率更高。齊義泉等(2003)[3]用WAVEWATCH-III模擬了1996年南海海浪場(chǎng),并對(duì)結(jié)果進(jìn)行了分析,他們采用的驅(qū)動(dòng)風(fēng)場(chǎng)是NCEP風(fēng)場(chǎng)。李明悝等(2005)[4]用QSCAT/NCEP混合風(fēng)場(chǎng)模擬了東中國海風(fēng)浪場(chǎng),參照日本氣象廳發(fā)布的浮筒數(shù)據(jù),證實(shí)風(fēng)場(chǎng)和計(jì)算所得有效波高與實(shí)測(cè)值符合得很好。
本文利用WAVEWATCH-III對(duì)具有復(fù)雜地形的近岸海域波浪演化進(jìn)行數(shù)值模擬。忽略風(fēng)場(chǎng)影響,采用江蘇科技大學(xué)所做的模型試驗(yàn)的地形與工況,與模型試驗(yàn)結(jié)果進(jìn)行分析對(duì)比,分析了誤差來源,驗(yàn)證了WAVEWATCH-III對(duì)近岸波浪場(chǎng)模擬的有效性,為應(yīng)用此模式進(jìn)行更深一步的波浪數(shù)值模擬提供了參考依據(jù)。
對(duì)于隨機(jī)海浪來說,海面的不規(guī)則變化可以用譜密度E描述,這種海浪譜通常稱作能量譜,可以表示為相參數(shù)的函數(shù)再考慮時(shí)間和空間的變化,就可以寫成其中為波數(shù)矢量,σ,ω分別為固有頻率和絕對(duì)頻率,分別代表地理空間和時(shí)間。參數(shù),σ,ω不是相互獨(dú)立的,它們通過波動(dòng)的頻散關(guān)系和多普勒方程建立聯(lián)系。WAVEWATCH-III模式選擇以波數(shù)矢量和方向θ為基本的參數(shù)組成譜但模式的輸出仍然采用以往模式的頻率和方向譜作為基本輸出,這兩種譜的轉(zhuǎn)換可以通過雅可比向前變換來實(shí)現(xiàn)。
假設(shè)水深和海流發(fā)生顯著變化的空間比波幅大,或者說,水深和海流的變化相對(duì)波浪來說是緩變,則準(zhǔn)非線性波浪理論成立,并存在下面的色散關(guān)系和多普勒效應(yīng)關(guān)系:
方程(1)-(3)即可解出波浪頻率參數(shù)隨時(shí)空的演化。
波譜E是所有獨(dú)立相位參數(shù)的函數(shù),同時(shí)隨著時(shí)間和空間發(fā)生變化;波譜也是方向的函數(shù)(方向譜),即也簡(jiǎn)單寫成和傳統(tǒng)的方向譜表達(dá)式一致。
不考慮海流因素時(shí),波包的能量是守恒的。Whitman等(1974)[5]曾提出,在背景流場(chǎng)不均勻時(shí),能量密度E是不守恒的,而波作用量密度N=E/σ是傳播守恒的。即
N代表波作用密度譜,S代表能量密度譜E中的所有源函數(shù),在WAVEWATCH-III中,
上式中各源項(xiàng)分別代表:線性輸入、風(fēng)—浪相互作用、耗散、波浪非線性相互作用、海底摩擦作用和碎波效應(yīng)(本計(jì)算不考慮風(fēng)場(chǎng))。
在直角坐標(biāo)系中,平衡方程(4)寫成:
其中:矢量Cg是群速度,由標(biāo)量群速度幅值Cg和方向θ決定,s是方向θ的一個(gè)坐標(biāo),m是垂直于s方向的一個(gè)坐標(biāo)。在笛卡爾坐標(biāo)系下(9)式是成立的。方程(7)-(9)中上標(biāo)一點(diǎn)代表對(duì)時(shí)間導(dǎo)數(shù),即x˙=算子是二維哈密頓算子。
方程(3)和(6)為主要的控制方程。求解該方程組就可以獲得時(shí)空中任一點(diǎn)的波譜。
地形采用模型試驗(yàn)對(duì)應(yīng)的實(shí)際地形。參照江蘇科技大學(xué)論文《近島礁生產(chǎn)生活平臺(tái)及浮式棧橋系統(tǒng)水動(dòng)力性能模型試驗(yàn)》[6],模型地形各測(cè)點(diǎn)位置說明見圖1所示,圖中w1~w18為監(jiān)測(cè)點(diǎn)位置。試驗(yàn)縮尺比1:36,量綱化的測(cè)點(diǎn)坐標(biāo)見表1。
將試驗(yàn)地形數(shù)據(jù)地形數(shù)據(jù)轉(zhuǎn)換成海底分布,同時(shí)將坐標(biāo)原點(diǎn)轉(zhuǎn)換到海平面的左下角。由于試驗(yàn)給出的地形分布網(wǎng)格間距較大且間距不等,需要進(jìn)行插值加密,得到了四種精度的插值方法(網(wǎng)格間距為4.5 m)得到的海底3D分布。由于原地形崎嶇不平,不同的差分精度得到的地形分布不同,總體上二階和四階得到的地形分布更不光滑,而一階精度較低,因此本文采用三階精度的插值結(jié)果作為地形分布。三階插值得到的海底3D分布如圖2。
圖1 w1~w18測(cè)點(diǎn)布置位置示意圖(單位:m) Fig.1 Arrangement positions of measuring points w1~w18(unit:m)
圖2 三階插值得到的海底3D分布Fig.2 Submarine 3D distribution obtained by different interpolation methods
表1 測(cè)點(diǎn)坐標(biāo)位置Tab.1 Arrangement positions of measuring points
圖3給出了地形網(wǎng)格間距對(duì)計(jì)算結(jié)果的影響。網(wǎng)格間距不同,插值得到的地形分布略有不同??梢钥闯鰀x=9 m的結(jié)果與dx=4.5 m的結(jié)果有一定的差距,而dx=4.5 m的結(jié)果與dx=2.25 m的結(jié)果接近。根據(jù)以往的仿真經(jīng)驗(yàn),波頻率離散101個(gè)點(diǎn)已經(jīng)足夠了,波方向離散72個(gè)點(diǎn)也已經(jīng)足夠了。本文計(jì)算的網(wǎng)格間距為dx=2.25 m根據(jù)以往的仿真經(jīng)驗(yàn),波頻率離散101個(gè)點(diǎn),波方向離散72個(gè)點(diǎn)。
圖3 不同網(wǎng)格距離對(duì)波高的影響Fig.3 The influence of different grid distance on the wave height
在模型初始化時(shí),WAVEWATCH-III模式提供了默認(rèn)設(shè)定與用戶設(shè)定兩種開關(guān),表2列出了本文所使用的模型設(shè)置參數(shù)。
頻率分段方法為:
其中分辨率Xσ取1.01,波數(shù)的空間分割數(shù)m取101。初始頻率σ0根據(jù)具體工況給出,方向分段數(shù)取72。
由于計(jì)算工況是不規(guī)則波,參照w18位置的試驗(yàn)波浪譜,計(jì)算出自定義譜,使得在波浪傳播入口w18點(diǎn)的波浪譜與試驗(yàn)波浪譜吻合。
自定義譜的主要輸入?yún)?shù)有:波高h(yuǎn)s和形狀函數(shù)其中ql角度逆時(shí)針增加,x方向?yàn)?,y方向?yàn)?0,而在WAVEWATCH-III中方向角度x方向?yàn)?70,y方向?yàn)?80。根據(jù)輸入?yún)?shù),WAVEWATCH-III計(jì)算出的能譜分布為:
其中:
表2 模型設(shè)置參數(shù)Tab.2 Model set parameters
在模型試驗(yàn)中,造波機(jī)的位置處于w17,造波機(jī)的實(shí)際造波能力與設(shè)計(jì)工況之間有一定誤差,所以在WAVEWATCH-III中,計(jì)算域入口處w17采用的工況為模型試驗(yàn)中w18處測(cè)得的工況經(jīng)過量綱變換得到。四種工況為不規(guī)則波工況,對(duì)于不規(guī)則波,入口處波譜根據(jù)試驗(yàn)波譜插值得到,頻率范圍選擇為0.05到0.25之間,其中橫坐標(biāo)為圓頻率(單位為rad/s),縱坐標(biāo)為能量(單位為cm2s/rad)。具體計(jì)算工況參數(shù)見表3。
表3 計(jì)算工況參數(shù)Tab.3 Calculated work condition parameters
本計(jì)算對(duì)檢測(cè)點(diǎn)進(jìn)行了監(jiān)控,并輸出計(jì)算定常后的監(jiān)控點(diǎn)的值,包括有義波高(m)、平均波長(zhǎng)(m)、平均周期(s)、平均方向角(°)、方向偏轉(zhuǎn)角(°)等。四種計(jì)算工況檢測(cè)點(diǎn)的信息以及與實(shí)驗(yàn)的比較(某些監(jiān)測(cè)點(diǎn)由于缺少試驗(yàn)數(shù)據(jù)而未給出),其中,波譜波高即為試驗(yàn)測(cè)得的的有義波高,波高絕對(duì)誤差等于波譜波高減去計(jì)算波高,波高相對(duì)誤差等于波高絕對(duì)誤差與計(jì)算波高之比,各個(gè)測(cè)點(diǎn)的計(jì)算結(jié)果在下表4-7中給出。
表4 工況EB401檢測(cè)點(diǎn)的信息以及與試驗(yàn)的比較Tab.4 Information of detection points and compared with the experiment both under work condition EB401
表5 工況EB402檢測(cè)點(diǎn)的信息以及與試驗(yàn)的比較Tab.5 Information of detection points and compared with the experiment both under work condition EB402
表6 工況EB403檢測(cè)點(diǎn)的信息以及與試驗(yàn)的比較Tab.6 Information of detection points and compared with the experiment both under work condition EB403
表7 工況EB501檢測(cè)點(diǎn)的信息以及與試驗(yàn)的比較Tab.7 Information of detection points and compared with the experiment both under work condition EB501
圖4-7給出了不規(guī)則波的四種工況的一維波譜對(duì)比,其中左邊為仿真結(jié)果,右邊為試驗(yàn)結(jié)果。
從圖4可以看出,仿真與試驗(yàn)在入口w18處的波譜完全一樣,在w2處的接近,在w6處的波峰值有一定的差別,w10處的波形也比較接近,但是試驗(yàn)中w14處的波譜除了在0.1 Hz處有一個(gè)峰值外,在0.2 Hz處又出現(xiàn)一個(gè)較大峰值,導(dǎo)致試驗(yàn)中w14處的波高偏高,為0.807 m,而仿真中的波高要低于入口w18處的波高,為0.526 m,因此導(dǎo)致w14處試驗(yàn)與仿真的波高誤差超過了50%。
圖4 工況EB401觀察點(diǎn)上的一維波譜對(duì)比Fig.4 Comparison of simulation and experiment of one-dimensional wave spectrum on observation points under work condition EB401
圖5 工況EB402觀察點(diǎn)上的一維波譜對(duì)比Fig.5 Comparison of simulation and experiment of one-dimensional wave spectrum on observation points under work condition EB402
從圖5可以看出,仿真與試驗(yàn)在w18處的波譜完全一樣,在w2處的接近,在w6處定性一致,而w10處和w14處的波譜在頻率大于0.2 Hz處仍較大。
從圖6可以看出,無論是仿真還是試驗(yàn),w2、w6、w10處的波譜與w18處的波譜接近。類似工況EB401,w10處試驗(yàn)波譜在頻率0.17 Hz處出現(xiàn)第二個(gè)較大峰值,導(dǎo)致試驗(yàn)與仿真的波高誤差超過50%。在w14處試驗(yàn)與仿真的波譜定性一致,但是試驗(yàn)中在頻率大于0.1 Hz范圍內(nèi)仍存在一定的能量,盡管試驗(yàn)和仿真在該處的波高接近,但是仿真波譜在0.08 Hz處的峰值要大于實(shí)驗(yàn)值。
從上面的分析可以看出,試驗(yàn)波譜與仿真波譜定性一致。但試驗(yàn)中w10和w14處的波譜存在兩個(gè)明顯的波峰,對(duì)于頻率較大的波峰仿真中未出現(xiàn)。由于第二個(gè)峰值的出現(xiàn),導(dǎo)致了試驗(yàn)與仿真的誤差較大。
從圖7可以看出,仿真中除了w14位置處的波譜明顯小于入口w18處的波譜外,w2、w6處的波譜與w18處很接近,特別是w2處的波譜與w18處幾乎重合。然而試驗(yàn)中w2處的波譜要明顯高于w18處的波譜,而w6處的波譜明顯小于w18處的波譜,因此導(dǎo)致w2處的波高實(shí)驗(yàn)與仿真有一定的誤差。可能認(rèn)為w6處的波譜不正確,試驗(yàn)報(bào)告中并未給出在該處測(cè)量的波高。w10和w14的波高對(duì)于試驗(yàn)和仿真是定性一致的,但試驗(yàn)中w10處的波譜在0.17 Hz處存在第二個(gè)峰值。
圖6 工況EB403觀察點(diǎn)上的一維波譜對(duì)比Fig.6 Comparison of simulation and experiment of one-dimensional wave spectrum on observation points under work condition EB403
圖7 工況EB501觀察點(diǎn)上的一維波譜對(duì)比Fig.7 Comparison of simulation and experiment of one-dimensional wave spectrum on observation points under work condition EB501
實(shí)際海域中的波為不規(guī)則波浪,本文對(duì)5種不規(guī)則波工況結(jié)果分析。從計(jì)算仿真結(jié)果來看,除了EB401的入射波波高較小時(shí)在w13和w14兩點(diǎn)上試驗(yàn)和仿真的誤差較大,其它工況和檢測(cè)點(diǎn)上的波高相對(duì)誤差絕大部分在10%左右。此外,當(dāng)入射波高較大時(shí),檢測(cè)點(diǎn)上波高的相對(duì)誤差較小。
在WAVEWATCH-III中,海底深度在以下5個(gè)方面對(duì)波浪演化存在影響:(1)波長(zhǎng)修正;(2)海底壁面摩擦;(3)非線性作用;(4)碎波;(5)限幅值。其中第一項(xiàng)是水波滿足的色散關(guān)系進(jìn)行波長(zhǎng)修正,第二至第四項(xiàng)通過模型源項(xiàng)來計(jì)算,第五項(xiàng)是水波假設(shè)理論的限制,在WW3計(jì)算過程中,有效波高限制在波長(zhǎng)的12%以下,而波長(zhǎng)又由色散關(guān)系與水深密切相關(guān),這5個(gè)方面的影響隨波高的增加而增加。而實(shí)際情況中,對(duì)于海底水深較淺而波高又較大時(shí),波浪往往會(huì)漫過淺攤,在淺灘處的水深會(huì)有一個(gè)明顯的抬升,形成壅水,而且此時(shí)波浪的傳播不一定是按照正弦波的形式傳播,水波理論假設(shè)不一定成立。因此采用WAVEWATCH-III仿真時(shí),對(duì)于海底水深較淺而波高又較大的情況存在一定誤差,而在實(shí)驗(yàn)與仿真的比較中,就會(huì)出現(xiàn)淺灘測(cè)點(diǎn)上的波高誤差偏大(如w13和w14點(diǎn))。在WAVEWATCH-III仿真中,如何合理對(duì)波長(zhǎng)進(jìn)行修正、如何定義淺灘的壁面摩擦系數(shù)、非線性模型和碎波模型、如何考慮淺灘中水深提高而改進(jìn)限幅值的方法需要進(jìn)一步的研究。
致謝
本研究由 973 項(xiàng)目 2013CB036101,國家自然科學(xué)基金(51379033;51379030;51239002;51221961),中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)資金(DUT15LK43,DUT15LK32),高科技船舶資助的項(xiàng)目(2016)22資助完成。
[1]Feng H,Vandemark D,Chapron B,et al.Use of a global wave model to correct altimeter sea level estimate[J].Geo Science and Remote Sensing Symposium,2004(4):2738-2741.
[2]Strauss D,Castelle B,Browne M,et al.Empirical estimation of near shore waves from a global deep-water wave model[J].IEEE Geo Science and Remote Sensing Letters,2006,3(4):462-466.
[3]齊義泉,朱伯承,施 平,等.WAVEWATCH-III模式模擬南海海浪場(chǎng)的結(jié)果分析[J].海洋學(xué)報(bào),2003,25(4):1-9.Qi Y Q,Zhu B C,Shi P,et al.Analysis of significant wave heights from WWATCH and TOPEX/Poseidon Altimetry[J].Acta Oceanologica Sinca,2003,25(4):1-9.
[4]李明悝,侯一筠.利用Quick SCAT/NCEP混合風(fēng)場(chǎng)及WAVEWATCH-III模擬東中國海風(fēng)浪場(chǎng)[J].海洋科學(xué),2005,29(6):9-12.Li M L,Hou Y J.Simulating wind-wave field of the East China Seas with QuikSCAT/NCEP blended wind and WAVEWATCH[J].Marine Sciences,2005,29(6):9-12.
[5]Whitman G B,Wiley J.Linear and nonlinear waves[M].New York:John Wiley&Sons,1974.
[6]凌宏杰,王志東.近島礁生產(chǎn)生活平臺(tái)及浮式棧橋系統(tǒng)水動(dòng)力性能模型試驗(yàn)[R].鎮(zhèn)江:江蘇科技大學(xué),2014.Lin H J,Wang Z D.Model experiment on hydrodynamic responses of floating platform near island and floating pier[R].Zhenjiang:Jiangsu University of Science and Technology,2014.
[7]Tolman H.User manual and system documentation of WAVEWATCH-III version1.18 NOAA/NCEP Tech[R].Note 166,1999:110.
[8]Tolman H,Booij N.Modeling wind waves using wavenumber-direction spectra and a variable wave number grid[J].Global Atmosphere and Ocean System,1998(6):295-309.
[9]Han S L.Evaluation of WAVEWATCH-III performance with wind input and dissipation source terms using wave buoy measurements for October 2006 along the east Korean coast in the East Sea[J].Ocean Engineering,2015(100):67-82.
[10]Tolman H.The numerical model WAVEWATCH:A third generation model for the Hind easting of wind wave son tides in shelf seas[R].Communication on Hydraulic and Geo technical Engineering,Delhi University of Technology,1989(2):89.
[11]Wittmann P A.Implementation of WAVEWATCH III at fleet numerical meteorology and oceanography center[J].OCEANS,MTS/IEEE Conference,2001(3):1474-1479.
[12]Hanson J L,Tracy B A,Tolman H L,Scott R D.Pacific hindcast performance of three numerical wave models[J].Journal of Atmospheric and Ocean Technology,2009(26):1614-1633.
[13]Tolman H,Balasubramaniyan B,Burroughs L D,et al.Development and implementation of wind-generated ocean surface wave models at NCEP[J].Weather and Forecasting,2002,17(4):311-333.
[14]梅強(qiáng)中.水波動(dòng)力學(xué)[M].北京:科學(xué)出版社,1983.Mei Q Z.Wave dynamics[M].Beijing:Science Press,1983.
[15]文圣常,余宙文.海浪理論與計(jì)算原理[M].北京:海洋出版社,1984.Wen S C,Yu Z W.Wave theory and calculation principles[M].Beijing:China Ocean Press,1984.
[16]高加云.波浪譜數(shù)學(xué)模型初步應(yīng)用研究[D].南京:河海大學(xué),2006.Gao J Y.Rudiment research in application of spectral models for wind-waves[D].Nanjing:HoHai University,2006.
[17]Miles J W.On the generation of surface waves by shear flows[J].Fluid Mech,1957(3):185-204.
[18]林建國,邱大洪,鄒志利.新型Buossinesq方程的進(jìn)一步改善[J].海洋學(xué)報(bào),1998,20(2):113-119.Lin J G,Qiu D H,Zou Z L.Further improvement of new Boussinesq-type equations[J].Acta Oceanologica Sinca,1998,20(2):113-119.
[19]郭衍游,侯一筠,楊永增,等.利用WAVEWATCH-III建立東中國海區(qū)域海浪同化系統(tǒng)[J].高技術(shù)通訊,2006,16(10):1092-106.Guo Y Y,Hou Y J,Yang Y Z,et al.To build a regional ocean wave data assimilation system of eastern China seas with WaveWatch III[J].Chinese High Technology Letters,2006,16(10):1092-1096.