岸礁地形上波浪增水和傳遞波的數值模擬研究

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

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

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

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

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

1 FUNWAVE-TVD數值模型簡介

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

1.1 控制方程與數值解法

FUNWAVE-TVD數值模型的控制方程為

η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和子網紊動混合項Rs。

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

1.2 邊界條件與造波

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

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

(6)

(7)

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

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

(8)

(9)

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

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

(10)

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

(11)

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

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

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

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

Rf=-Cd×uα|uα|

(12)

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

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

2.1 物理試驗簡介及數值模型設置

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

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

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

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

2.2 模型檢驗與參數敏感性分析

2.2.1 模擬結果對空間網格步長Δx的敏感性分析

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

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

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

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

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

2.2.2 模擬結果對底摩阻系數Cd的敏感性分析

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

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 模擬結果對破碎閾值H/d的敏感性分析

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

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

通過以上的岸礁波浪試驗模擬可以發現：FUNWAVE-TVD適用于岸礁地形上波浪傳播變形的數值模擬，能夠較合理地預測岸礁地形上波浪的時均水位變化、波高變化、礁后岸坡上波浪爬高等水動力過程。

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

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

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

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

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

4 結論

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

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

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

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