基于SWASH模型的珊瑚島礁波生環(huán)流模擬研究

陳松貴，王建平，高 駿，譚 召

(1.交通運輸部天津水運工程科學(xué)研究所，天津 300456；2.91053部隊，北京 100070；3.中船第九設(shè)計研究院工程有限公司，上海 200063；4.中交四航局港灣工程設(shè)計院有限公司，廣州 510220)

珊瑚礁是深海和淺海中均存在的一種珊瑚蟲骨骼沉積形成的巖體，我國廣闊的海域如東沙、西沙、中沙、南沙諸群島和臺灣島等區(qū)域擁有豐富的珊瑚礁資源。近年來，我國開展了一系列珊瑚島礁成陸工程，建設(shè)了諸如港口、機場、碼頭、燈塔、通信等相關(guān)設(shè)施，這些工程與設(shè)施在復(fù)雜的島礁海岸水動力環(huán)境下的穩(wěn)定性問題越來越受到學(xué)者們的關(guān)注[1]。不同于以往常見的近岸緩變地形的水動力環(huán)境，島礁附近水深劇烈變化且存在潟湖和裂口等。當(dāng)遠海波浪傳播到珊瑚礁時，由于礁前斜坡較陡、水深急劇變化，波浪在此發(fā)生強烈破碎，破碎波在礁坪上傳播消耗大量能量后繼續(xù)形成向前傳播的波浪，并且破碎波會在礁坪上產(chǎn)生部分波生流以回流的形式返回到外海[2]。

在這類環(huán)境下，研究波浪在珊瑚礁地形的傳播變形與增水等水動力特征，對工程建設(shè)與珊瑚礁保護具有重要意義。珊瑚礁地形上波浪傳播的研究方法涉及理論分析、現(xiàn)場觀察、物理模型實驗與數(shù)值模擬等。Gourlay[3]通過一系列水槽物理模型實驗研究了規(guī)則波作用下礁坪上的波生流問題，并給出了波生流經(jīng)驗計算關(guān)系式。Demirbile等[4]為了研究礁坪增水及岸灘波浪爬高隨不規(guī)則波及風(fēng)的關(guān)系，開展了風(fēng)浪水槽實驗研究，結(jié)果表明波浪爬高會因珊瑚礁而出現(xiàn)衰減。梅弢等[5]研究了50 a一遇和常年平均波高在概化珊瑚礁地形上的傳播過程，總結(jié)了礁坪位置外海波浪的傳播規(guī)律。Yao等[6]基于物理模型實驗研究了礁前斜坡坡度、礁坪水深與破碎波特征(破碎帶寬度、破波位置與破碎類型)之間的關(guān)系。柳淑學(xué)等[7]基于物理模型實驗研究了波浪破碎位置同波高變化的關(guān)系。姚宇等[8]在水槽實驗中證實了波浪增水、破碎受冠頂水深影響的現(xiàn)象。陳松貴、王建平等[9-12]對建有防浪堤的珊瑚礁陡變地形上波浪的傳播變形、越浪規(guī)律和波浪力特征進行了一系列研究，并總結(jié)了礁坪波浪增水出現(xiàn)的原因，即在礁緣處波浪破碎引起輻射應(yīng)力的急劇變化。Lowe等[13]討論了夏威夷Kaneohe灣珊瑚礁不同入射波況對波浪破碎和環(huán)流的影響，提出了一個理想的礁坪-潟湖-裂口環(huán)流模型。Yao等[14]在波浪水槽布置了帶裂口的珊瑚礁物理模型，模擬礁坪增水和波生流在裂口存在時的變化情況。近年來基于Navier-Stokes方程的非靜壓模型得到快速發(fā)展(SWASH、NHWAVE等)。此類模型省去傳統(tǒng)CFD模型模擬自由水面時采用的如VOF、level set等過程，直接假定自由水面是關(guān)于水平位置的單值函數(shù)[15]。Buckley等[16]分別應(yīng)用了XBeach模型、SWAN模形、SWASH模型模擬波浪在陡坡岸礁的傳播變形，結(jié)果表明在模擬礁坪上的波生流時非靜壓模型可以得到更為準確的結(jié)果。姚宇等[17]通過對珊瑚島礁附近的孤立波的傳播變形特征進行的數(shù)值模擬，證實了非靜壓模型NHWAVE相比于Boussinesq模型對垂向流速分布特征有著更高的模擬精度。鄒國良等[18]在SWASH模型中引入了渦粘模型以提高其模擬波浪破碎的能力，且使模型在較少的垂向分層下即可得到較好的結(jié)果。Zheng、Yao[19-20]等采用三維港池物理模型研究了礁坪-潟湖-裂口的環(huán)流特征和波浪分布。

本文采用SWASH模型，模擬了珊瑚島礁的三維港池實驗，并與物理模型試驗的波高和流速結(jié)果進行對比，驗證模型的準確性。最后給出了不同外海波浪條件下珊瑚島礁地形上波高、流速、增水和波生流的變化規(guī)律。

1 SWASH模型簡介

1.1 控制方程

SWASH(Simulating Waves till Shore)模型由荷蘭代爾夫特理工大學(xué)的Stelling、Zijlema和Smit等學(xué)者歷時近15 a開發(fā)完成[18-20]。模型考慮了波浪運動中的非靜壓影響，適用于深海至淺水條件下的波浪和水流運動的模擬，特別是對于復(fù)雜地形如三維區(qū)域下的波浪破碎所引起的增水和波生流有較好的適用性。SWASH的控制方程包括淺水方程、垂向動量方程和水平動量方程中附加的非靜壓項，在求解時把控制方程分解為靜壓項和非靜壓項兩部分分部求解[21]。其控制方程如下

(1)

(2)

(3)

(4)

式中：t為時間；x、y和z方向的速度分量分別為u、v、w；水平向粘滯系數(shù)為vH；垂向粘滯系數(shù)為vv；ξ為自由面；g為重力加速度；q為非靜水壓力。

1.2 邊界條件

自由面的壓力邊界條件為

q|z=ξ=0

(5)

當(dāng)垂向流速滿足式(5)的邊界條件時水平流速滿足

(6)

在SWASH中關(guān)注入流邊界位置，其邊界條件的設(shè)定，需要明確速度參數(shù)，同時考慮弱反射問題[22]

(7)

式中：ub為入流邊界的流速；ξb為邊界入射波的波幅；式(7)中邊界位置是通過正負號來描述的，如表現(xiàn)為+號，則入流邊界的位置處在南、西，如表現(xiàn)為-號，則位置處在東、北。SWASH中利用干濕處理技術(shù)來描述動邊界。在模型設(shè)置中給定最小水深值，在計算的單位時間步長中對網(wǎng)格點進行干濕判斷，當(dāng)計算的網(wǎng)格點水深低于最小水深時，將該網(wǎng)格點的水深設(shè)為最小水深，避免模型計算出現(xiàn)負水深引起的模型運行中斷。

1.3 數(shù)值求解

模型計算采用有限差分格式進行空間離散，垂向采用σ分層網(wǎng)格，變量通過交錯網(wǎng)格方式布置。非靜水壓力q定義在單元邊中心，如圖1所示。

圖1 網(wǎng)格垂向分層和變量布置示意圖

2 SWASH模型驗證

2.1 物理模型簡介

物理模型布置如圖2所示，區(qū)域長60 m、寬42 m，模型比尺為1:50。在距造波機34 m處設(shè)置坡度為1:8的斜坡模擬礁前斜面，斜面后接長度為14 m的水平平臺模擬礁坪，在礁坪中間存在寬度為6 m的裂口，礁坪后存在寬度為5 m的潟湖，最后在潟湖后設(shè)置坡度為1:3.3的礁后斜坡。實驗中采集波面使用的是47根LG1型電容式浪高儀(G1～G47)，由于該礁坪-潟湖-裂口系統(tǒng)物理模型關(guān)于裂口對稱，所以大部分儀器主要布置在裂口的一側(cè)。

圖2 模型布置圖(單位：m)

2.2 模型參數(shù)設(shè)置及驗證結(jié)果

本文主要研究規(guī)則波情況下珊瑚島礁地形的環(huán)流特性，對應(yīng)的驗證工況為：礁坪淹沒水深hr=0.04 m，外海入射波高H0=0.04 m，周期T0=2 s。波浪入射為垂直左邊界入射，為使模型計算穩(wěn)定在右邊界處設(shè)置5 m海綿層。垂向分層為4層，糙率設(shè)置曼寧粗糙系數(shù)n為0.01。

由圖3礁坪上測點的波面歷時過程對比可知，在礁緣波浪破碎處的數(shù)值模擬結(jié)果和物理模型結(jié)果吻合較好，說明SWASH模型對于波浪破碎時的波面變化可以較為準確的模擬。裂口和潟湖不同測點的波面歷時曲線數(shù)值模擬和物理模型趨勢基本一致，數(shù)據(jù)吻合效果較好。

圖3 波面數(shù)模結(jié)果和實測值對比

圖4分別對比了礁坪中心線(圖4-a)和裂口中心線處(圖4-b)在x方向的物理模型和數(shù)值模型的向岸平均流速，平均流速取造波穩(wěn)定后60個周期內(nèi)的平均?？梢钥闯鼋钙荷系南虬镀骄魉贁?shù)值模擬計算值與實測值吻合較好，在裂口中心線處數(shù)值模擬的最大向岸流速比物理模型的偏大，分析原因可能是物理模型是通過攝像機采集的表面浮子的運動軌跡來間接測量表面流場的流速，存在測量誤差，而數(shù)值模型是計算理想情況下的垂向平均流速，因此輸出的表面流速也會隨著垂向分層的數(shù)量變化有一定的影響，出現(xiàn)了數(shù)值模型偏大的結(jié)果。

4-a 礁坪中心線 4-b 裂口中心線

3 珊瑚島礁波生環(huán)流模擬

3.1 波浪增水分析

5-a 二維波高云圖 5-b三維波面分布圖

選取模擬結(jié)果中穩(wěn)定的波面過程求其平均值，得到平均水位(Mean water level, MWL)的空間分布，其相對于靜水位為正時則為增水，負值為減水。圖6是平均水位沿程變化情況，其中圖6-a為二維增水云圖，圖6-b為三維增水分布圖。由圖6可知在礁前斜坡上主要以減水為主，波浪在礁緣破碎后出現(xiàn)明顯的波浪增水并在靠近礁緣處達到最大，此后逐漸減小。礁坪上的增水向裂口處逐漸減小，這可能是由于在水平方向向裂口回流的影響所致。礁坪上的增水沿岸分布規(guī)律為向裂口處逐漸減小。

6-a 二維增水云圖 6-b 三維增水分布圖

圖7為不同外海波要素的礁坪向岸方向中線上沿程平均水位變化趨勢圖，可以看出礁坪上的沿程平均水位隨著波高的增大而增大。波高較大的工況中，在x=36.5 m附近達到最大減水，而入射波高較小的工況較晚達到最大減水(x=37 m)，原因是波高較小的工況在礁坪上的波浪破碎強度低于其他工況，且發(fā)生位置較晚。波浪破碎后沿程平均水位迅速增大，在x=40 m處達到最大增水，隨入射波高增長，最大增水遞增倍率在113%～154%。隨后沿著礁坪增水持續(xù)減小，同時在礁坪上產(chǎn)生向岸方向的壓力梯度，驅(qū)動向岸方向的波生流。由于波生流的動能在礁后斜坡轉(zhuǎn)化為壅高的勢能[9]，因此可以看到潟湖位置增水有增大的趨勢。礁坪上的沿程平均水位隨著周期的增大而增大。但是和改變波高相比，改變周期的平均水位的增長趨勢相對不明顯。礁坪上的沿程平均水位隨著礁坪水深的增大而減小，即在低水位時礁坪的淹沒深度要大于高水位，在工程建設(shè)中要考慮其影響。

7-a 規(guī)則波 T0=2 s，hr=0.04 m 7-b 規(guī)則波 H0=0.04 m，hr=0.04 m 7-c 規(guī)則波 H0=0.04 m，T0=2 s

3.2 環(huán)流流速分析

圖8 珊瑚島礁流場圖

由圖9可知，礁坪上流向裂口的沿岸流速隨著波高增大而增大，礁坪最大向岸流速隨著波高的增大而增大，但隨入射波高增大流速增長幅度減小。向岸流速和沿岸流速隨入射波周期T0的增大而增大，但影響較小。而流速和靜水深的關(guān)系較為復(fù)雜，存在某個靜水位(礁坪水深hr=0.04 m)使得流速達到最大值。

9-a 規(guī)則波 T0=2 s,hr=0.04 m 9-b 規(guī)則波 H0=0.04 m,hr=0.04 m 9-c 規(guī)則波 T0=2 s,H0=0.04 m

由圖10可知，潟湖處的向岸流速隨著波高的增大而增大，其最大值位置較為接近。和礁坪上流速的變化規(guī)律類似，潟湖的流速隨周期增大而增大，但波高影響較小。沿岸流和靜水深的關(guān)系較為復(fù)雜，在中等水位(礁坪淹沒水深hr=0.4 m)流速達到最大值，其原因可能與礁坪上類似，受水深變化引起的增水變化和平均阻力變化的相互作用的影響。

10-a 規(guī)則波 T0=2 s,hr=0.04 m 10-b 規(guī)則波 H0=0.04 m,hr=0.04 m 10-c 規(guī)則波 H0=0.04 m,T0=2 s

由于裂口中不同剖面上的流速差異較大，所以對比離岸流速最大的裂口中線剖面y=18 m處在不同工況中的流速情況。由圖11可知，離岸流速隨著波高增大而增大，且隨著靠近外海的口門，離岸流速增大，離岸流的流出距離增加。中等水位(礁坪淹沒水深hr=4 cm)時在裂口整個出流過程中達到最大流速值，與礁坪上流速隨水位變化規(guī)律相同。但不同水位時的流速變化過程有所差異，中水位時流速最大值出現(xiàn)在裂口偏后位置，其他水位最大值出現(xiàn)在礁緣附近。

11-a 規(guī)則波 T0=2 s，hr=0.04 m 11-b 規(guī)則波 H0=0.06 m，hr=0.04 m 11-c 規(guī)則波 H0=0.04 m，T0=2 s

4 結(jié)語

本文采用SWASH模擬了不同的外海波要素和水深對珊瑚島礁地形的波浪增水和波生環(huán)流特性及變化規(guī)律進行了分析。主要結(jié)論如下：

(1)波浪由深水傳播到礁前斜坡發(fā)生破碎，波高迅速衰減，礁坪邊緣靠近裂口的向岸剖面上沿程最大波高可達其他位置2倍以上。入射波高越大、靜水深越小的工況破碎后波高衰減更多，高水位工況波高破碎位置較晚，且破碎前波高增大比例更大。礁前斜坡發(fā)生波浪減水后行進約一倍深水波長達到最大增水，裂口的存在使得礁坪上靠近裂口位置的增水減小。礁坪上的波浪增水隨外海入射波高、周期的增大而增大，隨礁坪的淹沒水深增大而減小，與波高、水深相比，周期的影響較小。

(2)在增水正壓力和輻射應(yīng)力的驅(qū)動下礁坪上產(chǎn)生向岸流，在礁坪靠近裂口邊緣處和潟湖中發(fā)生方向旋轉(zhuǎn)，在裂口中形成離岸流通向外海，從而形成外海-礁坪-潟湖-裂口的水平環(huán)流系統(tǒng)。在固定裂口和潟湖寬度情況下，不同外海波浪入射要素和礁坪水深的工況的流速分析表明，礁坪、潟湖和裂口中流速隨入射波高增大而增大，而周期對礁坪上流速影響較大，即隨周期增大流速增大，對于澙湖和裂口中的流速影響較小。在不同礁坪水深中，受輻射應(yīng)力和底摩阻的相互作用影響，在中等水位時達到最大流速值。