孤立波對礁坪上直墻沖擊試驗和RANS數值模擬

肖 理，房克照，孫家文，2，范子初，劉忠波

(1. 大連理工大學 海岸和近海工程國家重點實驗室，DUT-UWA海洋工程聯合研究中心，遼寧 大連 116024； 2. 國家海洋環境監測中心，遼寧 大連 116023； 3. 國家海洋局秦皇島海洋環境監測中心站，河北 秦皇島 066002； 4. 大連海事大學 交通運輸工程學院，遼寧 大連 116026)

典型的珊瑚礁海岸由連接海床且較陡的礁前斜坡、延伸向岸灘的相對水平的礁坪組成，在礁前斜坡與礁坪相接處常存在隆起的礁冠[1]。島礁地形的存在，會引起水深的急劇變化，導致波浪在島礁上傳播時發生反射、折繞射、破碎和波能衰減等現象。這些變化后的波浪是近岸結構物的主要荷載，對水工建筑、軍事活動、港口和島礁建設都具有十分顯著的影響[2]。因此波浪在島礁地形上的傳播規律已成為海岸動力學研究的熱點。由于風暴潮、海嘯等海洋極端天氣引起的極端波浪事件是工程界最為重視的海洋水動力因素之一[1]。而島礁地形對波能衰減作用顯著，絕大部分的波浪在珊瑚礁上破碎并在礁坪上傳播時被耗散[3]，因此研究極端波浪下島礁附近波浪作用的水動力過程對于海洋防害減災具有重要的意義。由于在水動力方面存在相似性，孤立波一直被簡化用來研究海嘯波和島礁的相互作用問題[4]。

在物理試驗方面，Rober[5]早在2010年開展了孤立波在島礁地形上傳播的二維、三維水槽試驗，這些數據被廣泛用于數值模型驗證。Quiroga和Cheung[6]在試驗過程中引入底摩擦，研究了底床粗糙度對孤立波在島礁上傳播過程的影響。姚宇等[7]開展了一系列模型試驗，研究孤立波在不同島礁地形上的傳播變形和爬高過程，從入射波高、有無礁冠、潟湖長度和底床粗糙度等角度全方面展開了探討。這些研究基本集中在島礁剖面上水動力特征，對于孤立波對礁坪構上筑物沖擊的物理模型試驗較少。在數值模擬計算方面，通常采用淺水方程模型、非靜壓模型、Boussinesq模型以及CFD類模型進行[3]。Berger等[8]、Stefanakis等[9]、任智源等[10]利用非線性淺水方程研究了海嘯波與近岸島礁的相互作用，但非線性淺水方程基于靜壓假設推導，不能考慮波浪傳播演化過程中的色散效應。非靜壓模型考慮了動壓效應，可彌補這一不足。陳同慶和張慶河[11]利用改進后的SUNTANS三維非靜壓模型，對南海東北部海域內孤立波進行模擬；張其一等[12]運用SWASH 模擬了典型珊瑚礁地形上的波浪傳播過程，展示了非靜壓模型在珊瑚礁地形上良好的適用性。Boussinesq模型兼具色散性和非線性，可以較準確地模擬孤立波在島礁地形上的傳播過程。Roeber和Cheung[13]以及劉雨詩等[14]采用基于激波捕捉的Boussinesq方程成功地模擬了波浪在島礁上的傳播變形問題；孫志林等[15]利用Funwave-TVD模型準確地模擬了孤立波在岸礁地形的傳播和爬坡；房克照等[16]建立了具有間斷捕捉能力的 Boussinesq類波浪模型，并且用于模擬孤立波在島礁地形上的傳播，對比結果進一步表明了此類數值模型的精確度和適用性。但以上數值方法都不能有效模擬波浪破碎、相應的流場細節及波浪與構筑物相互作用，為克服這一缺陷，近年來一些學者通過采用CFD方法開展相關數值研究。高睿等[17]采用光滑粒子流體動力學(SPH)方法，研究了孤立波在斜坡地形上的傳播、變形和破碎過程。姚宇等[18]采用大渦模擬的方法(LES)成功地模擬了孤立波在珊瑚礁礁前斜坡上的淺化、礁冠附近的破碎以及破碎波在礁坪上的演化過程。

文中將利用水槽物理模型試驗和RANS數值模擬，研究典型島礁剖面上孤立波的傳播演化及其與礁坪后方直墻的相互作用過程。

1 試驗布置

1.1 水槽布置

試驗在大波流水槽中進行，水槽尺寸為69 m×2 m×1.8 m(長×寬×深)，圖1是試驗波浪水槽的布置圖。水槽側壁是玻璃，便于用相機記錄和直接觀察波浪傳播演變過程。用Goring[19]提出的造波方法，在水槽的一端用推板造波機產生孤立波。島礁剖面高0.4 m，前坡陡峭，坡度為1∶4，礁坪長6 m。坡面和礁坪表面由光滑的PVC板制成；在礁坪的末端安裝高0.6 m的直墻。礁體模型距造波機推板平衡位置38.89 m，由鋼架牢固支撐。礁冠也是由光滑的有機玻璃制成，可拆卸，高0.083 m，長1.55 m，兩側的坡度均為1∶4，安裝礁冠后形成較長的前礁坡，礁后斜坡形成潟湖。對于這兩種島礁地形，礁坪上最小水深為0.093 m，礁冠上最小水深為0.01 m。在試驗過程中，造波板處水深保持h=0.493 m恒定不變。以第四根浪高儀G4處測得的波高為入射波高(A=0.148 m)，無因次波高為ε=A/h=0.3，屬于非線性較強的孤立波。試驗重復三次，以獲得可靠的數據。

1.2 測量儀器

試驗一共采用33根浪高儀測量水槽沿程自由表面變化，它們的位置如圖1所示。G1-G3用于校核造波信號，G4-G7位于模型海床前端，用于監測入射波和反射波。G8恰好位于斜坡腳趾上方(即x=38.89 m)。G7-G21間隔0.2 m，覆蓋整個礁坡和部分礁坪，G21-G31間隔0.3 m，G32位于礁坪后方。最后一根浪高儀(G33)被放置在靠近垂直墻的位置，用來測量波浪爬高。四個ADV流速儀(聲學多普勒測速儀)分別安置在x=39.81、40.91、42.69、44.04 m處。ADV的探頭在礁面上方約0.04 m處。特別地，當地形安置礁冠時(x=40.49 m)，第二根速度傳感器ADV2被移動到x=42.07 m處。采用15個沿直墻分布的微型壓力傳感器記錄波浪對直墻的動態沖擊荷載，每個傳感器的直徑均為0.004 m，間距為0.02 m；精確地安置在直墻豎向軸線上(圖2是礁冠部分浪高儀布置圖和壓力測點分布圖)。所有的儀器在試驗前后都進行了校準，以確保線性和穩定性。自由表面、速度和壓力荷載的數據采集頻率分別為50 Hz、50 Hz和1 000 Hz。

圖1 試驗水槽和儀器布置示意Fig. 1 Layout of experimental flume and gague locations

圖2 礁冠處浪高儀和直墻上壓力傳感器布置示意Fig. 2 Layout of wave gauge over reef crown and pressure transducers on vertical wall

2 數值模型

采用OpenFOAM中的Waves2Foam求解器進行數值模擬。Waves2Foam是一個可以產生和吸收自由表面水波的工具箱。能夠產生孤立波、規則波、不規則波等各種波浪，該求解器的控制方程為連續方程和N-S方程，用VOF方法來追蹤自由表面，有限體積法離散控制方程。

2.1 湍流模型

采用RANS k-ε模型進行數值模擬，k-ε模型求解了兩個變量，k-湍流動能；ε-動能耗散率。OpenFOAM中標準的k-ε模型方程為：

(1)

(2)

k和ε的計算公式為：

(3)

其中，u是平均流速；I為湍流強度，I=0.16×Re-0.125；l為湍流尺度，一般取特征尺度的0.07倍，即l=0.07L；其余各常數取值分別為σk=1.00，σε=1.30，Cε1=1.44，Cε2=1.92，Cμ=0.09。

2.2 參數設置

理想化的數值波浪水槽應與物理水槽尺寸相同。但在實際模擬過程中，為節省計算資源和時間，下文只對物理島礁地形段進行了數值模擬，即從x=23.8 m到x=46.49 m，長度為22.6 m。由G4(x=29.13 m)記錄的波信號作為邊界入射波高。邊界條件設置如下：速度邊界條件在入口處設置成waveVelocity，頂部邊界設置為InletOutlet類型，出口設置成fixedValue，壓力在頂部設置為totalPressure，其余邊界處設置為zeroGradient；湍動能和耗散率在入口處設置為zeroGradient，底部和出口處fixedValue，由于建立的是一個二維數值水槽，前后兩個壁面條件設置為empty。

為確定網格尺寸，進行了敏感度分析。分別取方案一：Δx=0.02 m和Δy=0.008 m；方案二：Δx=0.01 m和Δy=0.004 m；方案三：Δx=0.005 m和Δy=0.002 5 m三組網格進行波面與網格尺度的敏感度分析，結果如圖3，方案三峰值明顯偏離試驗值，前兩組結果無明顯差異，考慮便于對波浪破碎的捕捉，最終確定網格尺度為Δx=0.01 m和Δy=0.004 m，網格在礁坪上沿x方向設置成漸變加密的形式，漸變系數為0.8；即最后一個網格的長度是第一個網格的0.8倍，以便更好地捕捉礁坪上的波浪破碎和直壁所受的動壓力。無礁冠網格示意圖見圖4，有礁冠時只在網格底部加入礁冠地形，無其他特殊處理。

圖3 波面變化與網格尺寸的敏感度分析Fig. 3 Sensitivity of wavefront change with mesh size

圖4 計算網格示意(礁坪上漸變加密)Fig. 4 Schematic diagram of computational grid (gradual refinement on the reef flat)

數值模型在40核的臺式工作站(2×20核Intel(R)Xeon(R)E5-2673 v4@2.30 GHz)上運行，內存32 GB RAM。模擬時間約為4小時，殘差精度控制在1×10-6階以內。各項數值模型參數設置如下：水的密度和運動黏度分別為ρ=1 000 kg/m3，ν=10-6m2/s；空氣的密度和運動黏度分別為ρ=1 kg/m3，ν=1.48×10-5m2/s。在礁冠頂端時雷諾數達最大(Re=dU/ν=11 000)。經計算給出了k-ε模型湍流動能k=2.94×10-6m2/s2和耗散率ε=1.27×10-7s-1。采用四邊形結構化網格對計算域進行離散。數值模擬的總時長為17 s，Courant數設置為0.45。時間步長在計算過程中自動調整，以滿足數值模擬過程的穩定性，其中最大時間步長設置為0.001 s。

3 數值結果分析

對比數值計算結果和試驗測量結果，包括自由面高程演變、流速、直墻上的波浪爬升和沖擊荷載。還分析了RANS k-ε模型計算得到的湍動能和渦量分布。數值計算結果和試驗測量數據在時空上都是同步的，以便進行比較，以孤立波波峰達到x=25.89 m的時刻定義數值計算時間零點(t=0 s)。

3.1 孤立波演化特征

3.1.1 無礁冠地形

試驗觀測到入射波在斜坡前保持完好的波形，在礁坪前端發生卷破，能量耗散，波幅減小，破碎波卷入空氣以湍涌的形式在礁坪上傳播并沖擊直墻，較大的流速使得直墻處產生極大的爬高，反射波浪二次破碎并以湍涌的形式向造波機方向傳播。

圖5(a)是特征時刻自由表面數值結果和試驗結果的對比圖。圖6(a)是對應時刻的湍動能和渦量場的數值模擬結果。由圖可見，孤立波傳播到礁前斜坡位置處時，由于淺水變形作用，波峰變陡并且略顯傾斜狀態(t=5 .5 s)，波形保持良好；在礁坪上水深更淺，波陡進一步增大，波浪明顯變形并開始發生卷破 (t=6.5 s)；產生的湍動能和渦量的量級大約是O(0.2 m2/s2)和O(15 s-1)。隨后波浪破碎，卷入大量空氣并以湍涌的形式向直墻傳播(t=7.6 s)。局部產生的湍動能和渦量在波頭處達最大值，并且通過渦旋向下后方進行傳播。最終波浪沖擊直墻并引起爬高(t=8.7 s)。最大湍動能和渦量在這一過程有所減小且分散分布，此時湍動能和渦量分布于大部分水底，之后在直墻前形成澭水，外海和礁坪上水位下降，直墻附近的水面有巨大抬升，形成反射波(t=9.3 s)，反射波二次破碎并形成湍涌向外海傳播(t=10.4 s，10.9 s)，這一過程湍動能和渦量聚集而增大。當波浪傳播到外海，礁坪上波面恢復平穩，水位下降，礁坪上產生持續的低頻波浪。整個過程數值計算結果和試驗結果吻合較好。

圖5 特征時刻自由表面的數值結果和試驗結果對比Fig. 5 Comparison of numerical results and experimental results of free surface

數值模擬觀察到8.7 s前后在水槽末端底部存在一個三角形區域，域內湍動能和渦量為零，這是由于破碎的波浪在水表面，且速度更快，湍動能通過渦旋向下和向后傳播，所以水中的湍動能和渦量分布要比水面的延遲出現(7.6 s)，另外在直墻和礁坪的阻隔下，使得水槽底部這部分水體所受擾動小，流速為零，形成靜水域，這也和試驗觀察相吻合。

3.1.2 有礁冠地形

圖5(b)是礁冠地形下特征時刻自由表面的數值結果和試驗結果。圖6(b)是對應時刻的湍動能和渦量場的數值模擬結果。礁冠延長了斜坡的長度，孤立波在礁前斜坡發生淺水變形的時間較無礁冠地形情況更長，且礁頂水深變淺。可以看到在礁冠頂端，波浪已經不具有穩定的波形，在礁冠末端卷破。此時湍動能幾乎為零，但是渦量聚集在一個點且量級達到了O(20 s-1)。可以認為處于臨界破碎狀態(t=6.0 s)；破碎點和無礁冠地形相比明顯提前。之后波浪傳播到潟湖，形成水躍，波浪破碎后卷入大量空氣并發展成湍涌在礁坪上傳播(t=6.7 s，7.0 s)；這一時刻波頭的湍動能最大值達到0.4 m2/s2，渦量最大值達到20 s-1。并通過渦旋向下向后傳播，從圖中可以看出渦量由聚集到分散的過程；隨后波浪沖擊直墻，湍動能和渦量急劇下降，分散于水底，在水槽末端存在一個靜水三角域。這和無礁冠地形情況類似(t=9.0 s)。t=10.5 s波浪發生反射并且反射波保持良好的波形，產生的湍動能達0.3 m2/s2，渦量只有5 s-1。對比無礁冠地形反射波破碎的情況，說明礁冠的存在對潟湖區極端波浪的產生有一定抵制作用。之后反射波向外海傳播，在礁頂和礁冠前坡均發生破碎(t=12.4 s，14.0 s)。

圖6 特征時刻湍動能和渦量數值模擬結果Fig. 6 Numerical simulation results of turbulent kinetic energy and vorticity

值得注意的是在t=7.4 s時，礁冠背面斜坡處出現了湍動能和渦量的極值。波形呈朝來波方向的卷舌狀，t=9 s波浪向來波方向破碎。從圖7可以觀察到，孤立波破碎后在礁坪上傳播過程中，大部分湍涌朝直墻傳播，與此同時礁冠背面產生一個反向的水流，部分波浪向來浪方向破碎并沖刷在礁冠背面，能量在礁冠背面耗散。這與姚宇等[20]試驗測得珊瑚礁礁冠會產生反向水流結果類似。

圖7 礁冠地形波浪傳播過程試驗照片Fig. 7 Snapshots of wave propagation process

3.2 波面和流速歷時曲線

圖8是無礁冠、有礁冠地形試驗和數值模擬波面歷時對比曲線。數值模擬能夠準確模擬孤立波主要演化過程。主峰的數值模擬和試驗結果取得較為精確的吻合，但在礁坪上，波形有破損(G26)。兩種地形上，浪高儀測點在采集到礁坪上波浪和反射波浪時出現了不同程度的高頻振蕩，然而在數值模擬中這些細節被光滑的曲線替代(G22、G26)，這是由于波浪破碎后形成的湍涌包含許多細小的渦旋，這些渦旋會造成電壓的擾動而在試驗采集中被記錄下來。而RANS模擬只能得到湍流平均物理量，捕捉不到這些小尺度渦旋。況且波浪破碎卷入空氣后是一個具有隨機性的過程，這對數值計算也是一個難點。另外，可以觀察到無礁冠、有礁冠地形入射波的尾波和反射波的尾波都不同程度的偏離試驗值，導致這一現象的可能原因是OpenFOAM數模造波和試驗室推板造波存在差異?？傮w而言，RANS k-ε封閉湍流模型能對孤立波傳播過程中自由表面進行較為準確的捕捉。

圖8 自由表面時間歷程對比曲線Fig. 8 Comparison of free surface time history

圖9比較了4個流速儀的水平速度分量的時間歷程。整個過程OpenFOAM都成功地捕捉到了波浪從外海傳播到礁坪和反射后的流速分布。其中無礁冠地形數值模擬結果吻合程度比礁冠地形高。ADV1位于斜坡中部，記錄礁前斜坡處波浪淺水變形過程中流速的變化，無礁冠地形下，波浪在坡前無反射，流速分布呈對稱。加上礁冠后，前礁坡有明顯的反射。流速下降到負值。ADV2位于礁坪前端，在無礁冠地形下，此時波浪變形但尚未破碎，流速分布不對稱，由于礁冠的存在，礁冠地形下ADV2移至x=42.07 m處，由圖5(b)和圖6(b)(t=7.4 s)以及圖7可知在礁冠背面產生反向的水流，與流速儀所測得的流速迅速降低到負值相吻合。而數模值沒有捕捉到這一變化的可能原因是空氣的混入影響到了數值模型的捕捉。ADV3以及ADV4數值模擬和試驗采集具有整體趨勢的相似性，但在細節上數值模擬并不能完全反映出來，這也可能是破碎波卷入空氣導致模擬不準確。從圖中可以看出，有礁冠地形下ADV3、ADV4所測的最大正向流速均小于無礁冠地形，這說明礁冠的存在使得礁坪上的流速減小?？赡艿脑蛴校航腹诘拇嬖谠龃罅斯铝⒉ǖ姆瓷洌瑫r改變了礁坪上的水深，且形成的反向水流在礁冠背面破碎，從而增強了孤立波破碎損耗的能量。

圖9 流速時間歷程曲線Fig. 9 Comparison of velocity time history

3.3 反射前沿程最大波面

圖10是無礁冠、有礁冠地形反射前沿程最大波面的試驗結果和數模結果對比圖。從圖中可以看出，無礁冠地形下波浪在礁前斜坡發生淺水變形，波面抬高卻不至于破碎，隨后最大波面有所降低，并在礁坪中部第二次抬高，波浪發生破碎，沿程最大波面變小；礁冠的存在使得波浪在礁冠處產生極大的抬高并破碎，之后整個過程沿程最大波面明顯減小。有礁冠地形的破碎點較之無礁冠地形提前，能量釋放更早，因此礁坪上傳播的最大波面要小于無礁冠地形的波面，且試驗記錄波浪在無/有礁冠兩種地形上反射時在直墻處的爬高分別為0.405 m 和0.347 m。礁冠的存在使得礁坪上的波浪變小，能大幅度減小波浪在直墻上的爬高，對潟湖、直墻起到一定的掩護作用。

圖10 沿程最大波面Fig. 10 Maximum wave height along wave flume

3.4 直墻上壓力分布

圖11是無礁冠、有礁冠地形直墻所受壓力的試驗和數模值對比結果。兩種地形下直墻所受沖擊過程相似，從P1到P11測點，壓力峰值先增大后衰減。整體而言，OpenFOAM對壓力的捕捉都是比較成功的，但并未準確地捕捉到P6～P8測點處的最大瞬擊荷載；從圖中可以看出，數值結果明顯小于試驗結果。最大誤差分別是45%(無礁冠地形P7測點)和54%(有礁冠地形P6測點)。分布在靜水面和水面以下的測點(P1～P5)，壓力先達到一個極值隨后又增大到最大值，然后衰減。位于水面以上(P6～P9)測點，壓力直接到達最大值，衰減后又增大到極值，之后衰減至零。這可能是波浪破碎沖擊直墻的位置處于P6～P9測點之間而產生的巨大瞬擊荷載，隨后產生的澭水導致了第二個極值的產生。圖12是所有測點最大壓力值，從下至上，最大壓力值有一個先增大后衰減的過程，結合G32測點測得的無/有礁冠地形湍涌最大高度為0.105/0.093 m，可知湍涌與直墻的作用位置分布在P5～P9區間。故在P6、P7、P8產生的瞬擊荷載很大,最大值出現在P6(礁冠)或P7(無礁冠)測點處，礁冠地形測得最大瞬擊荷載的測點位置要低一些，是因為礁冠的存在讓礁坪上最大波面低于無礁冠地形，沖擊的位置有所降低而導致的。對于試驗結果和數值結果都有個共同特點，除個別點外，有礁冠地形下的每個測點的壓力最大值都小于無礁冠地形下相應的測點。可以想象礁冠的存在，有助于減少直墻所受壓力。

圖11 壓力歷時曲線對比圖Fig. 11 Comparison of pressure time history

RANS k-ε模型能夠重現孤立波對直墻的沖擊過程。整體來看，數值模型對直墻測點的壓力變化有一個準確的描述，但當波浪破碎時，尤其是破碎和結構物共存時，在對瞬擊荷載的捕捉時，RANS k-ε模型模擬的結果明顯偏小。如圖12(b)、(c)是波浪沖擊直墻的瞬時照片和數值模擬的結果，可以看到沖擊過程中，波浪破碎嚴重且夾雜著大量的空氣，實際測量中瞬間的電壓變化可以被儀器記錄下來，但要在數值模擬中重現這一過程就顯得相當復雜且不準確。傳統的RANS模型只計算平均運動，而波浪破碎形成的細小渦旋且混入大量空氣，這是傳統的RANS模型無法處理的。對此，許多學者通過改進傳統的湍流模型來取得更高的精度，并結合試驗數據得到了很好的擬合效果，Brecht等[21]改進的k-ω SST湍流模型能夠準確捕捉孤立波的爬高；鄒學鋒等[22]結合傳統湍流模型，提出適用于珊瑚礁地形的k-ω SST湍流模型；Eltard和Fuhrman[23]發現傳統的湍流模型對湍動能預報偏高，改進并提出穩定多相流湍流模型(multiphase turbulence stabilization)。如果要得到更精準的數值模擬結果，可以提高網格質量并采用改進的湍流模型來進行數值模擬。

圖12 直墻所受最大壓力與試驗和數模結果照片Fig. 12 Maximum pressure on vertical wall and snapshots of experiment and simulation

4 結 語

采用物理模型試驗和數值模擬相結合的方式研究了孤立波在有礁冠和無礁冠的礁坪上的傳播過程及其對礁坪上直墻的沖擊作用。

對于文中的兩種島礁地形，發現入射孤立波都會在礁坪或礁冠上破碎，并發展成湍涌，在岸上傳播，最終沖擊直墻。結果表明，RANS k-ε模型能夠準確再現孤立波的主要演化過程以及波浪與直墻的相互作用，特別是能直接觀察到湍動能和渦量場的演變，這是試驗所不能得到的。不足之處是不能準確捕捉直墻最大瞬擊荷載，主要的一個影響因素是RANS不能捕捉破碎后的細小渦旋。且破碎的波浪卷入空氣使得整個沖擊過程具有很大的隨機性。從試驗數據和數模結果來看，由于礁冠的存在，波浪提前在礁頂破碎，礁坪上波高顯著降低，釋放能量較無礁冠地形更早；礁冠背面形成的反向水流對外海的來浪起到一定的阻擋作用，部分能量在礁冠背面釋放，使得礁坪上流速減小，從而減小了直墻上的爬升和最大動壓。所以礁冠的存在對近岸直墻、建筑物能起到一定的掩護作用。

當波浪破碎時，RANS k-ε模型模擬的結果顯得不夠精確，尤其是不能捕捉到瞬擊荷載，波浪破碎的研究一直是難點和熱點，傳統的湍流模型在處理波浪的破碎過程稍顯不足，為進一步研究孤立波與島礁地形上建筑物的作用，物理模型方面需要增加更多的試驗組次；數值模擬方面可以采用改進的湍流模型或者能夠捕捉小尺度渦旋的大渦模擬(LES)方法。