近岸波浪數(shù)值計算中兩種動邊界處理的比較

汪鴻，李訓(xùn)強，朱首賢，唐建，聶嶼

（1.解放軍理工大學(xué)氣象海洋學(xué)院 南京 211101；2.河海大學(xué)海洋學(xué)院 南京 210098；3.海軍南海艦隊 湛江 524000）

近岸波浪數(shù)值計算中兩種動邊界處理的比較

汪鴻1，李訓(xùn)強1，朱首賢2，唐建3，聶嶼1

（1.解放軍理工大學(xué)氣象海洋學(xué)院 南京 211101；2.河海大學(xué)海洋學(xué)院 南京 210098；3.海軍南海艦隊 湛江 524000）

在近岸波浪數(shù)值計算中，動邊界處理是不可避免的問題。文章基于Boussinesq方程的FUNWAVE-TVD模式，引入窄縫法對波浪水槽實驗進行數(shù)值模擬，比較窄縫法和干濕網(wǎng)格法的波浪數(shù)值計算結(jié)果；設(shè)計多種周期和波高的波浪數(shù)值試驗，分析兩種動邊界處理下數(shù)值計算結(jié)果的差異。結(jié)果表明，兩種動邊界處理對近岸波浪破碎的數(shù)值模擬有影響，對波浪漫灘區(qū)的計算影響尤為顯著。

動邊界處理；波浪破碎；波浪漫灘；海洋動力學(xué)；數(shù)值模擬

1　引言

波浪在由深水向淺水傳播過程中，由于水深變小，波陡不斷增大至極限，導(dǎo)致波浪發(fā)生倒卷和坍塌的現(xiàn)象，稱為波浪破碎［1-2］。波浪在破碎后向岸繼續(xù)傳播時，水陸分界線隨波峰向岸推進、隨波谷向?；赝?，在一個波動周期中灘地出現(xiàn)被海水淹沒和露出水面的過程，稱為波浪漫灘［3］。波浪漫灘的上界為波浪向岸爬高的最高位置，下界為波浪向?；赝说淖畹臀恢茫辖绾拖陆缰g的區(qū)域稱為波浪漫灘區(qū)。隨著科學(xué)技術(shù)的發(fā)展，人類對海洋資源的開發(fā)與利用顯著增加［4-5］，近岸海洋工程的建造與維護以及人員近岸活動需求也日益增多［6］，近岸海洋活動和工程建筑的安全保障成為較為突出的問題。近年的很多研究表明，波浪破碎是近岸區(qū)域波浪運動最為劇烈的物理過程，其對波浪預(yù)警預(yù)報以及近岸島礁的船只航行、工程建設(shè)和資源補給具有重要的安全威脅，因而一直受到廣泛關(guān)注［7-8］。近岸波浪破碎帶和波浪漫灘對近岸水動力有重要影響，是近岸泥沙運輸?shù)膬蓚€峰值區(qū)，是海灘泥沙侵蝕、搬運和地貌演變的重要機制［9-11］。一些研究表明，波浪漫灘區(qū)的泥沙輸運有時可占整個近岸區(qū)域沿岸方向泥沙輸運的50%［12-13］。深入研究近岸波浪演變有利于提高近岸波浪預(yù)警和保障水平，對維護近岸人員生命財產(chǎn)和國防安全具有重要意義。

數(shù)值計算是研究波浪近岸演變的重要手段，目前主要采用Boussinesq方程［14］和非線性淺水方程［15］進行近岸波浪的數(shù)值模擬研究。該類方程通過描述波動過程的水質(zhì)點運動來計算波浪，一般要求網(wǎng)格步長遠小于波長、時間步長遠小于波動周期，計算量龐大，適用于計算小范圍、短時間的波浪運動問題。這類模型在計算過程中，水陸間的動邊界處理是不可避免的問題，有些模型直接忽略近岸水線的移動，采用在一定水深處設(shè)置固壁的方法簡化處理邊界移動過程，而更多的波動質(zhì)點類模型采用的是干濕網(wǎng)格法［16-17］和窄縫法［18］處理波浪數(shù)值模擬過程中的動邊界問題。在干濕網(wǎng)格法計算時，隨著水陸交接位置的變化不斷更新計算邊界，并對有水區(qū)域重新進行網(wǎng)格劃分；該方法在計算中需要不斷重新劃分網(wǎng)格，機時耗費太大，在重新劃分網(wǎng)格時還容易出現(xiàn)病態(tài)的網(wǎng)格。Tao等［19］結(jié)合實際海岸中部分水體滲入灘地的情況提出窄縫法以處理近岸移動邊界，假設(shè)岸灘內(nèi)存在水流可以通過的窄縫，從而將動邊界區(qū)域轉(zhuǎn)化為始終有水的計算區(qū)域；該方法顯著提高計算效率，但是窄縫的存在造成波浪水體質(zhì)量和能量的損耗，這些效應(yīng)對波浪計算的影響一直備受關(guān)注，也有學(xué)者做了大量工作［18］。對兩種動邊界處理下近岸波浪數(shù)值計算的差異進行研究很有意義但比較少見，本文采用Boussinesq方程的FUNWAVE-TVD模式，分別引入窄縫法和干濕網(wǎng)格法，對近岸波浪數(shù)值計算結(jié)果進行比較和分析。

2　近岸波浪數(shù)值計算的動邊界處理

FUNWAVE-TVD模式是美國Delaware大學(xué)基于完全非線性Boussinesq方程建立的，該模式處理動邊界過程采用的是干濕網(wǎng)格法，對近岸附近的波浪傳播及演變具有較好的模擬效果［20］，但計算機時耗費較大。對于FUNWAVE-TVD模式的方程和處理有不少文獻介紹，本文不再贅述［14］。

窄縫法是目前處理動邊界的另一種主要方法，該方法把整個計算區(qū)域作為具有窄縫或可滲透的海床，把固體海底假設(shè)成具有窄縫或滲透邊界的活動型區(qū)域，使溢出的水位不會超過海灘高程（圖1）。Madsen等采用窄縫法對Carrier和Greenspan的理論模型進行數(shù)值模擬，發(fā)現(xiàn)其最大波高與理論解存在10%的誤差［21］。

圖1　窄縫法示意

引入窄縫法時，窄縫的控制參數(shù)為：

式中：κ為單位寬度上的窄縫寬；δ為最小縫寬；λ為窄縫的形狀參數(shù)；hz為窄縫起始的水深值；Λ為考慮窄縫影響后的等效水深；η為波面水位。本文引入窄縫法時，z*取為：

式中：zs為陸地高程。

3　近岸波浪數(shù)值計算結(jié)果比較

Hansen和Svendsen設(shè)計一組規(guī)則波在斜坡地形上傳播的水槽實驗［14］。水槽左邊的平底水深為0.36 m，右側(cè)斜坡坡度為Slope=0.029 2，波浪的入射波高H=4.3 cm、周期T=3.33 s。實驗主要測量波浪在斜坡上淺化傳播的波高變化（圖2）。

圖2　數(shù)值地形設(shè)置

本文采用FUNWAVE-TVD模式對Hansen和Svendsen的波浪水槽實驗進行數(shù)值模擬，數(shù)值模式的地形配置與波浪水槽實驗完全相同。造波源函數(shù)采用雙向造波法，造波源左端設(shè)置寬為1 m的海綿消波層以實現(xiàn)無反射造波，造波區(qū)位于斜坡左側(cè)10 m遠處。數(shù)值模型的網(wǎng)格距離為0.025 m× 0.10 m，時間步長為0.01 s。分別采用窄縫法和干濕網(wǎng)格法對波浪的動邊界進行處理。本文中窄縫形狀過渡參數(shù)為λ=20，窄縫寬度為δ=0.01 m，窄縫起始水深為hz=0.2 m，干濕網(wǎng)格判斷的臨界水深為0.001 m。模式計算20 s后，波面呈現(xiàn)很有規(guī)律的變化，波浪的波高基本保持穩(wěn)定。

將數(shù)值模式穩(wěn)定后30 s內(nèi)平均的波高模擬值與觀測數(shù)據(jù)比較，窄縫法和干濕網(wǎng)格法條件下的波高計算結(jié)果和實驗測量結(jié)果如圖3所示，其中橫軸為不同測點位置的水深，縱軸為各水深處波浪平均波高，“*”為數(shù)值計算的波浪破碎位置。結(jié)果表明，采用窄縫法和干濕網(wǎng)格法對FUNWAVE數(shù)值模型模擬波浪破碎的影響不大，即在兩種動邊界處理下，模型都能模擬出波浪由于地形淺化作用在近岸區(qū)域發(fā)生堆積并破碎的現(xiàn)象。波浪的波陡是波浪破碎的關(guān)鍵指標，在實驗測量數(shù)據(jù)中波浪破碎處波陡為0.029 6，采用干濕網(wǎng)格法和窄縫法計算的破碎波陡偏小，分別為0.027 3和0.024 5，證明計算的波浪破碎位置提前。在離岸較遠的區(qū)域，動邊界處理對波浪波高計算結(jié)果影響不大，隨著離岸距離的減小，波高的計算差異逐漸增大，在波浪破碎區(qū)和漫灘區(qū)波高模擬的差異尤其顯著。與干濕網(wǎng)格法計算結(jié)果比較，采用窄縫法處理動邊界下模擬的破碎位置前移，而破碎波高相對減小，計算的波浪在破碎區(qū)后波高也偏小，這種差異可能是由于窄縫的存在，在波浪向岸傳播過程中首先要填滿窄縫，從而引起水體質(zhì)量的損失并導(dǎo)致波浪能量的衰減，使波浪的破碎提前、破碎波高減小。與實驗數(shù)據(jù)比較，采用窄縫法和干濕網(wǎng)格法計算得到的最大波高分別減小0.82 cm和0.70 cm，破碎水深分別增大0.021 m和0.007 m；采用干濕網(wǎng)格法計算的破碎波高和破碎位置與實驗觀測結(jié)果更為符合，但波浪破碎后的計算波高顯著大于實測值。此外，與實驗數(shù)據(jù)相比，圖3中FUNWAVE模式波高的計算結(jié)果在破碎前偏大而破碎處偏小，Shi［14］和Kennedy［22］在采用Boussinesq方程計算波浪破碎時也出現(xiàn)破碎波高偏小的結(jié)果，Shi在模擬中提出該誤差是由于Boussinesq模型的網(wǎng)格數(shù)值耗散導(dǎo)致的［14］。

圖3　兩種動邊界下波高計算結(jié)果

波浪漫灘是近岸波浪運動的小尺度過程，水體的質(zhì)量和水質(zhì)點能量的變化對波浪漫灘過程的影響更為明顯。本文在FUNWAVE模式中采用窄縫法和干濕網(wǎng)格法處理Hansen和Svendsen實驗，分別計算波浪爬高的最高位置xh和回落的最低位置xl。采用窄縫法和干濕網(wǎng)格法處理動邊界，波浪最高和最低位置時水槽中波浪的瞬時水位如圖4所示?？芍捎谜p法計算的波浪爬高最高位置xh相對干濕網(wǎng)格法計算的xh偏低，而波浪回落最低位置xl則相差不大，很明顯是由于窄縫的存在，波浪在向岸爬升過程中水體的損失和能量的耗散導(dǎo)致的；圖4（a）中，由于窄縫的存在，在斜坡內(nèi)部仍然存在較小的水質(zhì)點波動，對波浪破碎后能量的損耗有明顯作用，這與實際觀測是較為相符的。

圖4　漫灘上下界位置兩種邊界計算水位分布

由于Hansen和Svendsen的實驗中沒有對波浪漫灘區(qū)進行觀測，本文采用目前工程應(yīng)用計算規(guī)則波爬高最多的Hunt經(jīng)驗公式計算Hansen實驗的爬高位置。Hunt公式為：

式中：R為波浪爬高的垂直距離；T和H分別為波浪入射周期和波高；β為坡度角。

計算Hansen實驗爬高位置所需的參數(shù)易由實驗條件得到。Hunt公式計算的爬高最大位置R= 0.237 m，轉(zhuǎn)化為沿波浪傳播方向上，波浪爬高區(qū)長為0.863 m。兩種動邊界處理得到的漫灘區(qū)范圍與Hunt經(jīng)驗公式結(jié)果的比較如表1所示，采用窄縫法和干濕網(wǎng)格法計算的波浪爬高區(qū)分別長0.675 m和1.30 m，與Hunt公式計算結(jié)果相比，波浪爬高區(qū)分別相差0.188 m和0.437 m，以Hunt公式計算結(jié)果為準的相對誤差分別為21.78%和50.64%，采用窄縫法計算誤差更??；結(jié)合圖3中波高的計算結(jié)果，這是干濕網(wǎng)格法在近岸漫灘區(qū)附近波高顯著偏大造成的。

表1　Hansen實驗漫灘區(qū)計算比較

4　不同入射條件下動邊界處理對波浪漫灘計算的影響

一般而言，波浪的入射波高和周期對波浪漫灘區(qū)的分布有較大影響。本文通過改變上述實驗的波高和周期，以Hunt公式計算的波浪爬高距離為標準，分別統(tǒng)計窄縫法和干濕網(wǎng)格法處理下不同入射條件得到的波浪爬高位置，并計算其與Hunt公式計算的爬高距離的相對誤差（表2和表3）。由于Hunt公式僅能計算波浪爬高位置，本文未討論兩種動邊界處理的波浪回落位置差異。

表2　入射周期不變（T=3.33 s）的波浪爬高

表3　入射波高不變（H=4.3 cm）的波浪爬高

表2保持波浪入射周期T=3.33 s不變，僅改變波浪入射波高；表3保持波浪入射波高H=4.3 cm不變，僅改變波浪周期。Hunt公式計算結(jié)果表明，隨著波浪波高和周期的增大，波浪的爬高距離隨之增加，采用窄縫法和干濕網(wǎng)格法都能較好地模擬出這種變化趨勢；由不同入射條件下兩種動邊界處理的爬高位置與Hunt公式計算結(jié)果的相對誤差可以發(fā)現(xiàn)，在本文選取的多種波浪入射條件下，采用窄縫法處理的爬高誤差顯著小于采用干濕網(wǎng)格法處理的爬高誤差，且窄縫法處理的相對誤差保持在20%左右，而干濕網(wǎng)格法處理的相對誤差變化較大。

5　結(jié)論

本文采用FUNWAVE-TVD模式對Hansen水槽實驗進行數(shù)值模擬，結(jié)果表明采用不同動邊界處理對近岸波浪破碎區(qū)的計算有較大差異。采用窄縫法計算的波浪破碎波高偏小且破碎位置相對前移，該法不僅能減少模型計算機時消耗，而且對波浪漫灘區(qū)計算的準確性也有一定的提高；對多種入射波高和周期條件下波浪爬高的計算結(jié)果進行比較，結(jié)果表明采用窄縫法處理與經(jīng)驗公式的相對誤差更小。在以后的工作中要加強對波浪爬高的現(xiàn)場觀測，通過數(shù)值模擬與實測資料的對比，進一步分析兩種動邊界處理對近岸波浪數(shù)值計算的影響。

［1］CHEN Qin，KIRBY J T，DALRYMPLE R A，et al.Boussineso modeling of wave transformation，breaking，and runup.II：2D［J］.Journal of Waterway Port Coastal& Ocean Engineering，2014，126（1）：39-47.

［2］QIN C，KIRBY J T，DALRYMPLE R A，et al.Boussinesq modeling of longshore currents［J］.Journal of Geophysical Research，2003，108（C11）：211-227.

［3］PULEO J A，BUTT T.The first international workshop on swash-zone processes［J］.Continental Shelf Research，2006，26（5）：556-560.

［4］鄭崇偉，潘靜，孫威，等.經(jīng)略21世紀海上絲路之海洋環(huán)境特征系列研究［J］.海洋開發(fā)與管理，2015，32（7）：4-9.

［5］鄭崇偉，李訓(xùn)強，高占勝，等.經(jīng)略21世紀海上絲路之海洋環(huán)境特征：風(fēng)候統(tǒng)計分析［J］.海洋開發(fā)與管理，2015，32（8）：4-11.

［6］鄭崇偉，李崇銀.中國南海島礁建設(shè)：風(fēng)力發(fā)電、海浪發(fā)電［J］.中國海洋大學(xué)學(xué)報：自然科學(xué)版，2015，46（9）：7-14.

［7］鄭崇偉，潘靜，黎鑫，等.1988-2009年中國海及周邊海域大浪頻率對ElNino的響應(yīng)［J］.海洋通報，2014（2）：140-147.

［8］ABADIE S C，GANDON S，GRILLI R，et al.3D numerical simulations of waves generated by subaerial mass failures：Application to La Palma case［C］//Proceedings of the 31st International Coastal Engineering Conference.Singapore：World Sci.，2008：1384-1395.

［9］HORN D P.Measurements and modelling of beach groundwater flow in the swash-zone：a review［J］.Continental Shelf Research，2006，26（5）：622-652.

［10］MASSELIND G，PULEO J A.Swash-zone morphodynamics［J］.Continental Shelf Research，2006，26（5）：661-680.

［11］BROCCHINI M.Integral swash-zone models［J］.Continental Shelf Research，2006，26（5）：653-660.

［12］WELLEN E V，BALDOCK T，Chadwick A，et al.STRAND：A Model for Longshore Sediment Transport in the Swash Zone［C］.Proceedings of the 27th International Conference on Coastal Engineering.Sydney：Coastal Engineering，2001：3139 -3150.

［13］Elfrink B，Baldock T.Hydrodynamics and sediment transport in the swash zone：a review and perspectives［J］.Coastal Engineering，2002，45（3）：149-167.

［14］SHI F，KIRBY J T，HARRIS J C，et al.A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation［J］.Ocean Modelling，2012，43-44（2）：36-51.

［15］BROCCHINI M，PEREGRINE D H.Integral flow properties of the swash zone and averaging［J］.Journal of Fluid Mechanics，1996，317（1）：241-273.

［16］張文靜，朱首賢，黃韋艮.衛(wèi)星遙感資料在湛江港風(fēng)暴潮漫灘計算中的應(yīng)用［J］.解放軍理工大學(xué)學(xué)報：自然科學(xué)版，2009，10（5）：501-506.

［17］呂新剛，喬方利，夏長水.膠州灣潮汐潮流動邊界數(shù)值模擬［J］.海洋學(xué)報，2008，30（4）：21-29.

［18］陶建華.波浪在岸灘上的爬高和破碎的數(shù)學(xué)模擬［J］.海洋學(xué)報，1984，6（5）：692-700.

［19］TAO Jianhua，LI Qingxue，F(xiàn)ALCONER R A，et al.Modelling and assessment ofwaterquality indicators in a semi-enclosed shallow bay［J］.Journal of Hydraulic Research，2001（6）：611-617.

［20］CHOI J，KIRBY J T，YOON S B.Boussinesq modeling of longshore currents in the Sandy Duck experiment under directional random wave conditions［J］.Coastal Engineering，2015，101：17-34.

［21］CARRIER G F，GREENSPAN H P.Water waves of finite amplitude on a sloping beach［J］.Journal of Fluid Mechanics，1958，4（1）：97-109.

［22］KENNEDY B，CHEN Q，KIRBY J T，et al.Boussinesq modeling of wave transformation，breaking，and runup I：1D［J］. Journal of Waterway Port Coastal &Ocean Engineering，2002，126（1）：39-47.

Comparison of Two Moving Shoreline Handling Method Results in near-Shore Wave Calculation

WANG Hong1，LI Xunqiang1，ZHU Shouxian2，TANG Jian3，NIE Yu1

（1.PLA University of Science and Technology，Nanjing 211101，China；2.Institute of marine，Hohai University，Nanjing 210098，China；3.Naval Nanhai fleets of PLA，Zhanjiang 524000，China）

Moving shoreline is an inevitable issue in numerical calculation nearing the surf zone. This paper laid slot method to the fully nonlinear Boussinesq model FUNWAVE-TVD to simulate flume experiment，then compared with the numerical result of slot and wet-dry method.Besides，series numerical calculations were designed to analyze the result of two shoreline handling method.The results showed that the impact of moving shoreline on wave shoaling and breaking calculation is obvious close，especially in the calculation of swash zone.

Moving-shoreline，Wave breaking，Swash，Ocean dynamics，Numerical simulation

P731

A

1005-9857（2016）09-0094-05

2016-04-11；

2016-07-25

國家自然科學(xué)基金（41206163，41376012，41076048）；中央高?；究蒲袠I(yè)務(wù)費項目（2011B05714，2014B06514）.

汪鴻，碩士研究生，研究方向為海洋動力學(xué)與數(shù)值模擬，電子信箱：wanghong_92@126.com

李訓(xùn)強，副教授，碩士，研究方向為海洋動力學(xué)與數(shù)值模擬，電子信箱：lixunqaing@sina.com