諸裕良，虞 琦，趙紅軍，宗劉俊

(河海大學(xué) 港口海岸與近海工程學(xué)院，南京 210098)

珊瑚礁是在熱帶和亞熱帶淺海、由造礁珊瑚骨架和生物碎屑組成的具有抗浪性能的海底隆起[1]，在中國臺灣島、海南島和南海諸島[2]也分布廣泛。根據(jù)礁體與岸線的關(guān)系可將其分為岸礁、堡礁和環(huán)礁[1]，其中岸礁沿岸生長并向海延伸。從地形上看，珊瑚礁海岸與常見的淤泥質(zhì)和沙質(zhì)海岸存在明顯的區(qū)別：淤泥質(zhì)和沙質(zhì)海岸的底坡一般較為平緩，波浪先在平緩的底坡上發(fā)生淺水變形，而后因水深的限制在淺水域破碎，破碎波在岸線附近形成沖泄區(qū)；而珊瑚礁由礁前帶、礁核帶和礁后帶組成[1]，礁前斜坡通常較陡，礁坪相對平坦，波浪在水深急劇變淺的礁前陡坡上發(fā)生淺水變形，隨后在礁緣附近破碎，破碎波在坦淺的礁坪上繼續(xù)傳播一段距離后生成新的行進波向礁后傳播。

國內(nèi)外學(xué)者對波浪在珊瑚礁地形上的傳播變形開展了大量的數(shù)值模擬研究，因良好的非線性性能Boussinesq方程數(shù)值模型在此方面得到了廣泛應(yīng)用。Skotner和Apelt[3]利用擴展型Boussinesq方程[4]，研究了規(guī)則波在潛礁上傳播過程中時均水位的變化，指出Boussinesq方程可較準確地模擬波浪破碎引起的增減水，僅在強非線性情況下計算的增水值偏小。Nwogu和Demirbilek[5]采用任意水深層Boussinesq方程[6]模擬了岸礁地形上不規(guī)則波的傳播變形，研究了波高和時均水位的沿程變化、以及由于非線性波波相互作用和波浪破碎等引起的波譜變化。Yao等[7]基于完全非線性Boussinesq方程數(shù)值模型[8]，模擬了不同岸礁地形上的波浪傳播變形，討論了礁前的坡度和剖面形狀對時均水位和波高的影響。Roeber和Cheung[9]采用有限體積法，建立了擴展型Boussinesq方程和淺水方程聯(lián)合計算模型，模擬了孤立波在岸礁上的傳播變形過程，利用試驗和觀測數(shù)據(jù)證實了模型能夠再現(xiàn)岸礁地形上波浪破碎前后的波面時間過程和譜形變化情況。房克照等[10]建立了高階Boussinesq 方程和淺水方程聯(lián)合計算模型，模擬了孤立波在潛礁上的傳播和破碎過程。

以往的研究多關(guān)注入射波浪條件[11,13]、礁坪水深[12-13]、礁前斜坡[7,12]、礁冠[13]等因素對珊瑚礁上波浪傳播變形的影響。但隨著珊瑚礁海岸的開發(fā)，礁頂工程建設(shè)日益增多，這改變了礁坪的長度和礁后岸坡的坡度，諸如此的地形變化對波浪增水和傳遞波的影響尚鮮有報道。為此，本文采用完全非線性Boussinesq數(shù)值模型FUNWAVE-TVD[14]，對岸礁地形上波浪的傳播變形進行研究，重點在于測試FUNWAVE-TVD對岸礁地形上波浪傳播變形的計算性能，研究礁坪長度和礁后岸坡坡度對波浪增水和傳遞波的影響。

1 FUNWAVE-TVD數(shù)值模型簡介

Kirby等[15]采用有限差分法建立了完全非線性Boussinesq方程數(shù)值模型FUNWAVE。隨后Shi和Kirby等[14]在此基礎(chǔ)上對控制方程和數(shù)值解法等做了若干改進，得到了FUNWAVE-TVD數(shù)值模型。

1.1 控制方程與數(shù)值解法

FUNWAVE-TVD數(shù)值模型的控制方程為

ηt+▽·M=0

(1)

uα,t+(uα·▽)uα+g▽η+V1+V2+V3+R=0

(2)

(3)

A=▽·(huα)

(4)

B=▽·uα

(5)

式中：V1、V2、V3是Boussinesq方程的色散項，R是擴散項和耗散項，包括底摩阻項Rf和子網(wǎng)紊動混合項Rs。

FUNWAVE-TVD數(shù)值模型采用有限體積和有限差分相結(jié)合的空間離散方法，其中，色散項使用常規(guī)的有限差分法處理，其他項使用有限體積法求解。模型采用高階MUSCL-TVD格式處理通量項和一階導(dǎo)數(shù)項，結(jié)合HLL近似黎曼求解器計算數(shù)值通量，從而使模型具備激波捕捉功能；采用具有強穩(wěn)定性的、步長三階Runge-Kutta法進行時間積分，使用CFL準則限制時間網(wǎng)格步長以保證計算收斂。

1.2 邊界條件與造波

FUNWAVE-TVD可以模擬固壁邊界、吸收邊界和動邊界。其中，固壁邊界采用比較常用的處理方法：令固壁處的法向速度、波面的法向梯度和切向速度的法向梯度都為零，即v=0，ηy=0，uy=0。

吸收邊界通過設(shè)置海綿層來實現(xiàn)，數(shù)值處理方法為在x、y方向的動量方程中分別加入人工阻尼項

(6)

(7)

式中：w1、w2、w3分別對應(yīng)牛頓冷卻項、粘性阻尼項和海綿層過濾項。海綿層寬度通常取2～3個波長。

動邊界采用干濕網(wǎng)格法處理，令干網(wǎng)格邊界面上的法向通量為零，同時將鏡面邊界條件運用于四階MUSCL-TVD格式和干網(wǎng)格色散項的離散，并對干網(wǎng)格處的波速進行修正

(8)

(9)

模型入射波邊界使用Wei等[16]提出的源函數(shù)造波法，在廣義形式的線性化Boussinesq方程和淺水方程中加入造波源函數(shù)項f(x,y,t)后，連續(xù)方程和動量方程可分別表達為

ηt+h▽u+α1h3▽2(▽u)=f(x,y,t)

(10)

ut+g▽η+αh2▽2ut=-g▽P

(11)

式中：α1和α為方程參數(shù)，取不同的值可得到不同形式的Boussinesq方程；P為壓力分布。

1.3 波浪破碎與底摩阻的處理

FUNWAVE-TVD通過關(guān)閉Boussinesq方程中的非線性項和色散項，將Boussinesq方程退化為非線性淺水方程來模擬波浪破碎。模型選擇波高水深比H/d作為方程的切換變量，對于任一網(wǎng)格單元，當(dāng)H/d≥0.8[17]時，則認為波浪發(fā)生破碎，此時非線性項和色散項不參與計算。

模型考慮底摩阻的方法是在動量方程中加入底摩阻衰減項

Rf=-Cd×uα|uα|

(12)

式中：Cd為底摩阻系數(shù)，Shi等[14]認為該參數(shù)在實驗室地形下不敏感，可取為0或者0.001。

2 岸礁地形上波浪傳播變形的數(shù)值模擬

2.1 物理試驗簡介及數(shù)值模型設(shè)置

為研究風(fēng)浪對礁坪波浪增水和岸坡波浪爬高的影響，Demirbilek等[18]在密歇根大學(xué)的風(fēng)浪水槽中開展了80多組次的波浪水槽模型試驗，包括波浪試驗、風(fēng)試驗和風(fēng)浪試驗。本文將利用其中的部分試驗，對FUNWAVE-TVD關(guān)于岸礁地形上波浪傳播變形的計算效果進行檢驗。

Demirbilek等[18]的試驗設(shè)置如圖1：試驗斷面的礁前斜坡為三段式復(fù)合斜坡，礁坪長4.8 m，高0.5 m，礁后岸坡坡度為1:12。沿程布置了九根電容式波高儀(G1～G9)，其中G1～G3位于坡腳前，用以分離入反射波，G4～G6、G7～G9分別位于礁前斜坡和礁坪，用以記錄沿礁波面時程的變化，礁后岸坡上布置了一根1 m長的電容式爬高儀測量波浪爬高。試驗采用不規(guī)則入射波，入射波譜型為JONSWAP譜，譜峰升高因子γ=3.3。波高儀和爬高儀的采樣頻率fs均為20 Hz，造波開始后立即采樣，采樣時間為15 min。

圖1 試驗設(shè)置圖(m)Fig.1 Experimental setup(m)

表1 入射波浪條件表Tab.1 Incident wave conditions

2.2 模型檢驗與參數(shù)敏感性分析

2.2.1 模擬結(jié)果對空間網(wǎng)格步長Δx的敏感性分析

空間網(wǎng)格步長Δx會影響模型的數(shù)值耗散和計算效率。為了測試模型計算結(jié)果對Δx的敏感性，針對Test33(礁坪無水)和Test45(礁坪淹沒)兩組試驗，分別以Δx為1/16Lp、1/32Lp和1/64Lp進行了計算，圖2給出了不同Δx時有效波高和時均水位的沿程變化，為方便比較，試驗值也繪于圖中。結(jié)果顯示：數(shù)值結(jié)果對Δx比較敏感，Δx較大時沿程波高計算值較試驗值明顯偏小，這是模型數(shù)值耗散過大所致，特別是在礁前斜坡段，波高迅速減小，數(shù)值耗散導(dǎo)致的波高衰減顯著強于淺水變形引起的波高增大，此時模型無法合理描述淺水變形和波浪破碎引起的波高變化。從時均水位看，Δx較大時減水現(xiàn)象不明顯且礁坪增水幅度明顯偏小，根據(jù)Tait[19]的研究，波浪的增減水與破碎波高成正相關(guān)，因為Δx較大時的數(shù)值耗散使破碎波高的計算值偏小，所以減水值與增水值亦偏小。

2-a1 Test33模擬與試驗有效波高沿程變化圖2-a2 Test45模擬與試驗有效波高沿程變化圖

2-b1 Test33模擬與試驗時均水位沿程變化圖2-b2 Test45模擬與試驗時均水位沿程變化圖

2-c1 Test33地形、水位示意圖2-c2 Test45地形、水位示意圖圖2 Test33和Test45模擬與試驗有效波高、時均水位沿程變化圖Fig.2 Predicted and measured significant wave height and mean water level variation for Test33 and Test45

為進一步探究Δx對波浪破碎模擬效果的影響，統(tǒng)計了不同Δx時的波浪破碎位置，并用相應(yīng)灰色豎線標記在圖2-a中。據(jù)此可知，在Δx較大時，模型預(yù)測的破碎位置后移且破碎范圍縮窄，這是因為模型以波高水深比作為破碎閾值，過大的數(shù)值耗散使波高計算值偏小，所以波浪在水深更淺時才發(fā)生破碎，從而造成破碎位置后移。圖3和圖4分別給出了波浪破碎前四根波高儀處模擬波譜與試驗波譜的比較，再次證明Δx過大帶來的數(shù)值耗散會使波能偏小，同時從圖中可以進一步看出數(shù)值耗散主要發(fā)生在不規(guī)則波中頻率大于譜峰頻率的高頻段。

圖3 Test33模擬與試驗波譜圖Fig.3 Predicted and measured wave spectra at selected gauges for Test33

比較圖2中有效波高和時均水位的計算結(jié)果和試驗值可知，在Δx取1/64Lp時，二者吻合良好，數(shù)值結(jié)果能較好地反映波浪在岸礁上傳播過程中波高和時均水位的變化。并且從圖3～圖4中也可以看出，在Δx取1/64Lp時，模型能較好地描述不規(guī)則波中各頻率的波能分布，在二倍頻(兩倍的入射波譜峰頻率)處模擬值與試驗值也基本吻合，說明模型具有較高的計算精度。

圖4 Test45模擬與試驗波譜圖Fig.4 Predicted and measured wave spectra at selected gauges for Test45

2.2.2 模擬結(jié)果對底摩阻系數(shù)Cd的敏感性分析

底摩阻系數(shù)Cd的取值影響數(shù)值計算中的底摩阻耗散量，Shi等[14]認為在實驗室條件下，模型對Cd不敏感，建議取0或0.001。但Nwogu和Demirbilek[5]在運用物理試驗對Boussinesq模型進行檢驗時發(fā)現(xiàn)：Cd的取值會影響礁坪段波浪模擬的精度。為確定適合本次試驗計算的Cd值并測試模擬結(jié)果對Cd的敏感性，對Test33和Test45兩組試驗再次進行計算。Demirbilek試驗為了減小摩阻對試驗結(jié)果的影響，使用了相對光滑的聚氯乙烯板制作地形，但實際操作時水槽內(nèi)無法做到絕對光滑，因此本文模型計算測試了Cd為0.001和0.006兩種情況。圖5對有效波高和時均水位的計算值與試驗值進行了比較。從圖中可以看出，模擬結(jié)果對Cd很不敏感，Cd取0.001和0.006對礁前斜坡上的波高和沿程時均水位幾乎沒有影響，僅小幅度地影響了礁坪上的傳遞波高，且Cd越大，傳遞波高越小。對比計算值與試驗值可知，Cd取0.001時，數(shù)值結(jié)果與試驗數(shù)據(jù)吻合良好，模型能更好地模擬岸礁地形上波高和時均水位的沿程變化。

5-a1 Test33模擬與試驗有效波高沿程變化圖5-a2 Test45模擬與試驗有效波高沿程變化圖

5-b1 Test33模擬與試驗時均水位沿程變化圖5-b2 Test45模擬與試驗時均水位沿程變化圖

5-c1 Test33地形、水位示意圖5-c2 Test45地形、水位示意圖圖5 Test33和Test45模擬與試驗有效波高、時均水位沿程變化圖Fig.5 Predicted and measured significant wave height and mean water level variation for Test33 and Test45

2.2.3 模擬結(jié)果對破碎閾值H/d的敏感性分析

FUNWAVE-TVD以波高水深比H/d作為破碎模型啟停的判據(jù)變量，Tonelli和Petti[17]建議將H/d的閥值取為0.8。為驗證該值在岸礁地形上的適用性以及測試模擬結(jié)果對H/d閥值的敏感性，對表1中的八組試驗進行了模擬。諸裕良等[20]通過波浪水槽試驗分析得出：珊瑚礁地形上波浪破碎時H/d大致分布在0.7～1.1之間，故計算中H/d的閥值分別設(shè)為0.6、0.8、1.0和1.2。因為波浪破碎后會在礁坪產(chǎn)生增水，破碎波在礁坪上傳播一段距離后生成新的行進波，隨后行進波在礁后岸坡上爬落，所以H/d的取值對波浪破碎以及破碎后的一系列水動力過程都會產(chǎn)生影響，因此以礁坪上的傳遞波高、增水和礁后岸坡上波浪爬高的相對誤差作為H/d閥值的評價標準。具體統(tǒng)計方法為：計算各組次下波高儀G8和G9處(G8和G9位于礁坪中后段，此時波浪已基本完成破碎，傳遞波高和增水都達到了相對穩(wěn)定的狀態(tài))的有效波高和增水，算出模擬值與試驗值的相對誤差，將波高和增水在兩個波高儀處的相對誤差分別平均，得到各組次的傳遞波高相對誤差和增水相對誤差。同時統(tǒng)計各組次下礁后岸坡上波浪爬高R2%的相對誤差。圖6以破波相似參數(shù)ξ0(ξ0=tanα/(H0/Lp0)0.5)，深水波高H0和深水譜峰波長Lp0均基于微幅波理論求得，H0=Hs/ks，ks為淺水變形系數(shù)，Lp0=gTp2/2π，tanα取礁前斜坡的平均坡度)為橫坐標，給出了H/d各閾值下傳遞波高、增水和爬高的相對誤差以及三者整體的相對誤差圖。從圖6-a～6-c可以看出，傳遞波高和增水的相對誤差都在40%以內(nèi)，爬高誤差稍大，但除了一組誤差較大外，其余各組也都在65%以內(nèi)。從圖6-d中可以看出模擬結(jié)果對H/d較為敏感，H/d為0.8時，數(shù)值計算結(jié)果與試驗值的整體誤差最小，所以模擬岸礁波浪傳播變形時H/d取0.8數(shù)值計算結(jié)果更準確。

6-a 傳遞波高相對誤差統(tǒng)計圖6-b 增水相對誤差統(tǒng)計圖6-c 爬高相對誤差統(tǒng)計圖6-d傳遞波高、增水和爬高的總相對誤差統(tǒng)計圖圖6 傳遞波高、增水、爬高以及整體相對誤差統(tǒng)計圖Fig.6 The statistics of the relative error of transmitted wave height, wave setup, wave runup and overall

通過以上的岸礁波浪試驗?zāi)M可以發(fā)現(xiàn)：FUNWAVE-TVD適用于岸礁地形上波浪傳播變形的數(shù)值模擬，能夠較合理地預(yù)測岸礁地形上波浪的時均水位變化、波高變化、礁后岸坡上波浪爬高等水動力過程。

3 礁坪長度和礁后岸坡坡度對波浪增水和傳遞波的影響研究

隨著珊瑚礁海岸的開發(fā)和資源利用，礁頂工程建設(shè)日益增多，這改變了礁坪的長度和礁后岸坡的坡度，為掌握此類變化對礁坪增水和傳遞波的影響，利用FUNWAVE-TVD模型開展數(shù)值模擬研究。以圖1的岸礁地形為原型，每次模擬僅改變礁坪長度B或礁后岸坡坡度m。為考慮極端情況，選取Demirbilek等[18]的岸礁波浪試驗中入射波高較大的組次Test31(有效波高Hs=9.2 cm，Tp=2 s，hr=1.6 cm)作為入射波浪條件。

3.1 礁坪長度B對波浪增水和傳遞波的影響

為消除水深對模擬結(jié)果的影響，基于微幅波理論將譜峰波長換算至深水，設(shè)置礁坪長度B分別為Lp0、2.5Lp0、5Lp0、7.5Lp0和10Lp0。不同礁坪長度時的有效波高和時均水位的沿程變化如圖7所示，可見：隨著礁坪長度減小，礁坪增水變化很小，但傳遞波高會逐漸增大。為更為直觀地考察增水和傳遞波與礁坪長度的關(guān)系，圖8給出了波浪增水、傳遞波高和爬高隨B的變化圖。為減少礁后岸坡上波浪淺水變形及破碎對傳遞波和礁坪增水的影響，以礁后岸坡坡腳前1 m處的有效波高和時均水位作為傳遞波高值和礁坪增水值，以礁后岸坡上累計頻率為2%的爬高R2%作為岸坡爬高值，由圖8可知：當(dāng)B從10Lp0減小到5Lp0，爬高幾乎不變，傳遞波高和增水會略有波動，但整體變化較小，其中波高增大了10%左右，增水減小了2%左右；當(dāng)B繼續(xù)減小時，傳遞波高和爬高不斷增大且增大速度越來越快，增水略有減小，B從5Lp0減小到Lp0過程中傳遞波高增大了20%，爬高增大超過70%，增水僅減小了6%?？偟膩砜?，礁坪長度縮短對增水幾乎沒有影響，但會使傳遞波高和爬高增大。

7-a 不同礁坪長度下有效波高沿程變化圖7-b 不同礁坪長度下時均水位沿程變化圖7-c 不同礁坪長度下地形、水位示意圖圖7 不同礁坪長度下有效波高、時均水位沿程變化圖Fig.7 Variations of significant wave height and mean water level over reef profile with different reef flat lengths

8-a 不同礁坪長度下傳遞波高變化圖8-b 不同礁坪長度下增水變化圖8-c 不同礁坪長度下爬高變化圖圖8 不同礁坪長度下傳遞波高、增水和爬高變化圖Fig.8 Variations of transmitted wave height, wave setup and wave runup with different reef flat lengths

9-a 不同礁后岸坡坡度下有效波高沿程變化圖9-b 不同礁后岸坡坡度下時均水位沿程變化圖9-c 不同礁后岸坡坡度下地形、水位示意圖圖9 不同礁后岸坡坡度下有效波高、時均水位沿程變化圖Fig.9 Variations of significant wave height and mean water level over reef profile with different beach slopes

10-a 不同礁后岸坡坡度下傳遞波高變化圖10-b 不同礁后岸坡坡度下增水變化圖10-c 不同礁后岸坡坡度下爬高變化圖圖10 不同礁后岸坡坡度下傳遞波高、增水和爬高變化圖Fig.10 Variations of transmitted wave height, wave setup and wave runup with different reef beach slopes

3.2 礁后岸坡坡度m對波浪增水和傳遞波的影響

根據(jù)人造岸坡和天然岸坡的坡度范圍，設(shè)置礁后岸坡坡度m分別為1:1.5，1:3，1:6，1:12，1:25，1:50。不同礁坪長度時的有效波高和時均水位沿程變化如圖9所示：隨著礁后岸坡坡度增大，礁坪增水變化很小，但傳遞波高會略有增大。圖10給出了礁后岸坡坡腳前1 m處的傳遞波高和增水值隨m的變化、以及岸坡上的爬高R2%隨m的變化圖。據(jù)此可知：當(dāng)m小于1:12時，傳遞波高和爬高隨著m的增大而增大，m從1:50增加到1:12的過程中，傳遞波高增大了18%，爬高增大了19%，增水增大不到1%；隨著m繼續(xù)增大，傳遞波高和增水都趨于不變，爬高先急劇減小后下降的速度變小，陳國平等[21]關(guān)于不規(guī)則波爬高的模型試驗研究也得到了這一規(guī)律，即爬高隨坡度先增大后減小。由此可見，礁后岸坡坡度增大不影響礁坪增水，在坡度較緩的范圍內(nèi)，坡度增大會使傳遞波高和爬高增大但增大的有限，隨著坡度進一步增大，傳遞波高不再受影響，但爬高會先顯著減小后減小的趨勢變緩。

4 結(jié)論

本文利用完全非線性Boussinesq方程數(shù)值模型FUNWAVE-TVD，對珊瑚岸礁地形上的波浪增水和傳遞波進行了數(shù)值模擬研究，得到以下結(jié)論：

(1)模擬岸礁波浪試驗的數(shù)值結(jié)果與試驗值吻合良好，說明模型可以對波浪在岸礁上的傳播變形過程進行有效計算，適用于岸礁地形上波浪傳播變形的數(shù)值模擬研究。

(2)數(shù)值模型的敏感性計算分析顯示：在利用該模型模擬岸礁波浪傳播變形時，建議空間網(wǎng)格步長取1/64Lp，破碎閾值取0.8，實驗室條件下的底摩阻系數(shù)取0.001。

(3)岸礁礁坪長度和礁后岸坡坡度改變基本不影響礁坪增水；當(dāng)礁坪長度小于5Lp0后，繼續(xù)減小礁坪長度會使傳遞波高和爬高增大且增大速度不斷加快；礁后岸坡坡度在坡度較緩的范圍內(nèi)增大時，傳遞波高和爬高隨之增大，隨著坡度進一步增大，傳遞波高不再受影響，但爬高會先急劇減小后減小速度變緩。