島礁地形上斜坡堤前破碎波流場數值模擬研究

黃晶鑫，楊 氾，潘軍寧，王登婷，劉清君

（南京水利科學研究院 河流海岸研究所，江蘇 南京 210024）

島礁是一種岸坡非常陡峭的海岸地形，波浪從深水傳播至島礁時，水深發生急劇變化，使波浪發生淺水變形并破碎，隨后在平坦的礁坪上繼續向前傳播。波浪破碎過程中產生的波動流場，對島礁底部沖淤變化、島礁周圍生物環境及防波堤護底和護面結構等產生很大影響。因此研究島礁地形上波浪破碎過程中的波動流速分布具有重要意義。

國內外學者采用現場觀測、物理模型試驗及數值模擬等手段對波浪在島礁上的傳播變形進行了大量研究。何棟彬等［1］采用非線性Boussinesq 型方程波浪模型，進行了規則波和隨機波在礁坪上傳播的數值模擬。譚安平等［2］應用FLOW-3D 數值軟件開展了珊瑚礁地形下防波堤波浪爬高特性研究。諸裕良等［3］采用非線性Boussinesq 型方程數值軟件FUNWAVE-TVD 模擬了礁坪上波浪破碎過程與增水特征。劉林平等［4］基于FUNWAVE-TVD 對孤立波在理想化三維島礁地形上的傳播及爬坡開展了現場尺度的平面二維數值模擬，結果表明：珊瑚礁的存在總體上可有效降低島嶼四周孤立波的最大爬坡高度；入射波高、礁坪水深、礁坪寬度、珊瑚礁糙率是影響珊瑚島礁四周孤立波爬坡分布的主要因素。

除上述基于歐拉法的數值模型外，拉格朗日型無網格粒子法——光滑粒子流體動力學法（smoothed particle hydrodynamics，簡稱SPH）也常被用于波浪破碎的模擬。該方法于1977 年分別由Lucy［5］，Gingold 和Monaghan［6］獨立提出，最初應用于解決天體物理學中的星系碰撞、流星沖擊等問題，之后發展至自由表面流［7］、多相流［8］、流固耦合［9］等問題的數值求解當中。溫鴻杰［10］基于SPH方法建立了一個能夠精確描述OWC式波能轉換裝置附近水動力特性的湍流模型，并改進了控制方程中的耗散項，克服了在進行長時間的數值模擬時自由液面處流體粒子向上漂移的現象。任冰等［11］建立了基于SPH 法的二維數值波浪水槽，模擬了規則波對透空式結構物的沖擊作用，并通過黎曼解和CSPM 相結合的方法對連續方程和動量方程進行了修正。溫鴻杰等［12］基于SPH 模型采用周期性邊界較好地模擬了規則波在深海島礁上的傳播變形過程以及礁坪上波浪增水空間分布。Dalrymple和Rogers［13］模擬了波浪越浪和淺灘破碎等問題并驗證了SPH方法的適用性。Altomare 等［14］模擬了采用復雜護面塊體的防波堤在波浪作用下的波浪爬高，驗證了SPH 方法在模擬該類問題中的可行性。Khayyer 等［15］在不可壓SPH 法的基礎上進行了修正，改善了波浪破碎的模擬效果。Lowe等［16］研究了島礁上不規則破碎波的流場。但目前尚缺乏對島礁上破碎波作用下斜坡堤前波動流場的相關研究。

鑒于波動流場尤其是近底流速對堤前沖刷預測和護底設計的重要性，基于SPH 模型建立數值波浪水槽，開展規則波作用下斜坡堤對島礁地形上波浪傳播變形和破碎波流場影響的數值模擬研究。首先模擬無防波堤時波浪在礁坪上的破碎過程，將礁緣前后不同位置處波動流速最大值垂線分布與波浪水槽物理模型試驗測量結果進行比較，以驗證SPH 模型的準確性和適用性；之后在數值水槽中添加斜坡堤，模擬分析建堤后礁坪上波動流場及波高的變化特征。

1 數值模擬方法

數值模擬采用基于SPH 方法的開源程序DualSPHysics［17-18］，該程序應用了CUDA 技術以支持GPU 運算，具有較高的性能和較好的適用性。

1.1 控制方程及其離散格式

弱可壓SPH法中的控制方程為N-S方程，連續性方程和動量方程為：

式中：ρ為流體密度；t為時間；u為速度矢量；P為應力張量；g為重力加速度；Γ為耗散項。

在SPH 法中，一個粒子在某一時間點的物理量通過使用核函數被離散化為該粒子周圍支持域內臨近粒子的物理量的求和，物理量的梯度采用核函數的梯度來表達。

核函數需要滿足一定條件［19］，主要包括：1）緊支性，在支持域外時，函數值恒為0；2）局域性，在支持域內時，函數值恒大于等于0；3）衰減性，粒子距離增大時，函數單調遞減。采用的核函數為Wendland［20］提出的五階樣條函數：

式中：αD為歸一化常數，在二維問題中αD= 3 (2πh2)，在三維問題中αD= 15 (16πh2)；q=r/h，無量綱量，為初始粒子間距與光滑長度的比值。

利用式（3）和式（4）將式（1）、式（2）中的物理量離散化后即可得到粒子離散形式的控制方程。在弱可壓SPH法中，粒子的質量為恒定值，計算過程中粒子密度發生變化［14］。

δ-SPH法在連續性方程中引入了人工黏性項，提高了計算的穩定性。采用了δ-SPH法的連續性方程為：

式中：α為人工黏度系數，研究［14，18］表明α= 0.01 為保證計算穩定性及避免非物理性振蕩的最小值；-cij=0.5(ci+cj)為平均聲速；μij=huij·rij(rij2+η2)為黏性系數，其中η2= 0.01h2；uij=ui-uj、rij=ri-rj分別為粒子i與粒子j處的速度矢量差和位置矢量差。

在研究礁坪上波浪破碎問題時，人工黏性理論中采用的人工黏度系數會造成過大的數值耗散［21］，文中計算根據邊界層黏性理論（laminar viscosity）［22］及亞粒子尺度湍流模型（sub-particle scale turbulence）［13］對其進行了優化。采用了邊界層黏性理論和亞粒子尺度法后的動量方程離散形式：

式中：ν0為邊界層中的運動黏度；τi、τj分別為i、j處的應力張量。

當采用弱可壓假定時，粒子處壓力依據粒子處密度通過狀態方程計算得出，數值計算中采用的時間步長由流體中的聲速計算確定［23］。為了使計算時間步長保持在合理范圍內，弱可壓假定中調整了流體的可壓縮性，從而可以采用較低的流體聲速。為了保證計算穩定性，密度變化需限制在1%以內，因此聲速需要大于10倍的最大流體速度。狀態方程為［24］：

式中：ri為粒子i的位移量。

1.2 邊界條件

水槽中的固邊界采用動態邊界條件法（dynamic boundary condition，簡稱DBC）［25］。該方法將結構邊界離散為粒子后，邊界粒子視為同樣滿足流體粒子的控制方程，并參與流體粒子控制方程的計算，但邊界粒子的位置不發生改變。

造波方法采用推板式造波板法，基于二階規則波造波理論，造波板位移量為：

式中：e(t)為造波板位移量；S0=H m1為造波板沖程；ω= 2πT為波浪圓頻率；?為初始相位；H為波高；k= 2πL為波數；d為靜水水深。

在使用推板式造波板法的過程中應用主動消波法可以消除二次反射波的影響。主動消波法（active wave absorption system，簡稱AWAS）［26］通過檢測造波板前的實際波面升高與理論值的差異來調整造波板的運動，從而對波面進行修正。造波板運動修正方程為：

式中：ηR(t)為理論波面升高和實際波面升高的差值；ηI(t)為理論波面升高；ηSPH(t)為SPH 計算所得的造波板前的實際波面升高；UR(t)為造波板速度的修正量。

水槽末端消波采用阻尼消波方式。阻尼消波通過修改粒子行進過程中的速度實現消波：

式中：v為修改后的粒子速度；v0為修改前的粒子速度；Δt為時間步長；β為可調節的速度削減參數；x為當前粒子位置；x0、x1分別為消波區的起始位置和結束位置。

1.3 時間積分

時間積分方法采用對偶位置Verlet 法（symplectic position verlet scheme），該方法具有二階精度，并且在采用了人工黏性項時仍具有時間可逆性和對稱性。當采用顯式時間積分法時，時間步長取決于庫朗數（courant-friedrichs-lewy，簡稱CFL）、壓力項和黏性擴散項。

1.4 模型設置

采用上述SPH 方法建立島礁地形上破碎波作用的二維數值水槽模型，圖1 所示為未設防波堤的數值水槽布置示意圖。水槽整體尺寸為50 m×1.5 m（長×高）；礁前靜水水深0.6 m，礁坪靜水水深0.1 m；坡前水平段長5.0 m，斜坡坡度為1∶1，長0.5 m，高0.5 m。水槽最左側為推波板，產生向右行進的規則波。試驗中波高H為0.085 m，周期T為1 s。水槽末端為人工阻尼消波區，以消除水槽右邊界反射波的影響。

圖1 未設防波堤的數值水槽布置示意Fig.1 Schematic diagram of a numerical wave flume in the absence of a slope breakwater

1.5 模型驗證

為驗證SPH 模型在模擬島礁地形上的波浪傳播變形方面的準確性，在波浪水槽中進行了物理模型試驗，開展對比研究。圖2 所示為物理模型試驗水槽布置示意圖。水槽尺寸為60.0 m×1.8 m×1.6 m（長×寬×高）。水槽一端配有推板式造波機，在兩端均配有消浪緩坡及消浪裝置。礁緣距造波機26 m，礁坪高0.5 m，試驗水深及波要素與數值模型相同。水槽中選取3條流速測線，分別位于島礁斜坡中間處、礁緣前、礁坪上，1#、2#和3#測線距礁緣距離分別為-0.25，-0.10，0.80 m，各測點垂向高度見表1。另在礁緣前0.50 m、礁緣處及礁緣后0.50 m 處分別設置A、B、C 三個波高儀。流速測量采用Vectrino 小威龍聲學多普勒點式流速儀（ADV），波高測量采用LG1型電容式波高儀，采用DS30型64通道浪高水位儀系統進行數據采集。

表1 模型驗證測點位置Tab.1 The location of measurement points for model validation 單位：m

圖2 物理模型試驗水槽布置示意Fig.2 Schematic diagram of a physical wave flume

將數值模擬得出的各測點正向及反向水平波動流速最大值的垂線分布與物理模型試驗結果進行對比（見圖3），可見數值模擬的波動流速與試驗測量結果總體符合較好，3#測線處數模計算的流速幅值略小。由于波浪破碎導致能量損耗，礁坪上3#測線處波動流速顯著小于礁緣前1#和2#測線處波動流速。圖4給出了1#、2#和3#測線近底水平波動流速變化過程驗證情況，圖5 為3 個波高測點處的波面過程對比圖，從圖中可以看出，SPH 模型的近底流速和波面過程模擬結果與物理模型試驗結果呈現良好的一致性。數值模擬與物理模型試驗在A波高儀處的波高分別為8.72、8.89 cm；在B波高儀處的波高分別為6.12、6.31 cm；在C波高儀處的波高分別為4.65、4.92 cm。

圖3 數值模擬結果與物理模型試驗結果中各點水平波動流速正、反向幅值對比Fig.3 Comparison of the forward and reverse amplitudes of the horizontal fluctuation velocity of the numerical simulation results and the physical model test results

圖4 不同測點處近底水平波動流速變化過程對比Fig.4 Comparison of near-bottom horizontal fluctuation velocity variation at various measuring points

圖5 不同測點處波面升高過程對比Fig.5 Comparison of wave surface elevation processes at various measuring points

由圖3～5 可見：在礁坪前，A 點處波面過程已表現出波峰前波面上升速度大于波谷前波面下降速度、波峰陡峭、波谷平緩及波谷歷時大于波峰歷時等非線性前傾波特征，但1#和2#測線近底水平波動流速仍為線性正弦波形態，正、反向水平波動流速最大值垂線分布均表現為自上而下逐漸減小的趨勢。在礁坪上，由于波浪破碎，波面和近底水平波動流速過程的非線性特征更為顯著：B 和C 點處波峰高度顯著大于波谷深度，波峰歷時顯著小于波谷歷時；3#測線近底正向流速最大值小于反向流速最大值，正向流速歷時也小于反向流速歷時，正向水平波動流速最大值由上向下逐漸減小，反向水平波動流速最大值則逐漸增大。

2 礁坪上斜坡堤對波動流場的影響

為研究斜坡堤對礁坪上波動流場的影響，分別采用SPH 數值水槽進行了有、無斜坡堤情況下的數模計算。計算中采用規則波，波浪要素、礁前水深等參數均與前文相同。圖6所示為設有防波堤的數值水槽布置示意圖。斜坡式防波堤設置在礁緣后0.8 m處，其坡度為1∶2，長6 m，高3 m，其余尺寸如圖6所示。

圖6 設有防波堤的數值水槽布置示意Fig.6 Schematic diagram of a numerical wave flume with a slope breakwater

圖7為有堤和無堤情況下不同時刻島礁波動流場對比圖。由圖7可見：有斜坡堤時，波浪在防波堤前發生反射，波動水體在堤坡上呈現周期性爬高、回落現象，波浪反射及爬高后回落水流對防波堤前流場造成顯著影響。t=0時，前波在破碎后波峰傳播至堤腳處，后波波峰在礁緣前發生劇烈變形，前緣變得陡峭，有堤與無堤時礁緣前至礁坪中段的波浪形態基本相同，有堤時出現波浪沿堤坡上爬現象。t=0.25T時，有堤時前波在堤坡上進一步爬高，后波波峰在礁坪上發生卷破，破碎點波高大于無堤時。t=0.50T時，后波在破碎后形成破后波，有堤時前波爬高水體回落，后波受反射波影響破碎紊動更為強烈，在波面下形成漩渦。當t=0.75T時，后波波峰傳播至堤前，有堤時波峰較無堤時顯著抬高，并再次發生卷破。總體上看，防波堤對礁坪上的波動流場形態產生很大影響，使波浪破碎更加猛烈，產生漩渦和紊動等復雜流體運動現象。

圖7 有堤和無堤時的島礁部分流場示意Fig.7 Schematic diagram of part of the reef flow field with and without slope breakwater

為具體研究斜坡堤對島礁流場的影響，對有、無堤時礁坪上不同位置處最大水平波動流速進行對比分析。從礁緣（圖6中5.5 m 處）至堤腳每隔0.1 m 設置一條垂線，每條垂線布置6個計算點，計算點垂向高度分別為0.560、0.540、0.530、0.520、0.510、0.505 m。圖8 為有、無堤時各點水平波動流速正、反向幅值對比圖，從圖中可見流速變化隨水平位置呈現顯著波動：有堤時，距礁緣0.1～0.3 m 及0.7 m 處正向流速較無堤時減小，其他位置正向流速增大，其中距礁緣0.2 m 處最大減小50%左右，距礁緣0.5 m 處約增大100%；有堤時反向流速除距礁緣0.1～0.2 m 及0.6～0.7 m 處略有增大或與無堤時接近外，其余點均顯著增大，其中距礁緣0.4 m處最大增大約90%。

圖8 有堤和無堤時各點水平波動流速正、反向幅值對比Fig.8 Comparison of forward and reverse amplitudes of horizontal fluctuation velocity with and without slope breakwater

圖9 為有堤和無堤時波高及近底水平波動流速幅值沿程分布圖。無堤時，由于波浪破碎后水體強烈紊動導致能量損耗，波高和波動流速均從礁緣處向礁緣后方逐漸減小。有堤時，由于堤前波浪反射的影響，入射波和反射波疊加導致波高和近底流速幅值均呈現沿程起伏波動現象：波高在距堤腳0.649 倍波長和0.108倍波長附近出現峰值，流速減小；波高在距堤腳0.455 倍波長附近出現谷值，流速則增大。總體上符合波浪部分反射時堤前波高和流速變化規律，但在堤前半倍波長范圍內（5.85～6.30 m區間段）有堤時反向近底流速峰值增幅較大，顯著大于波高峰值增幅，且正向和反向流速峰值出現位置略有不同。

圖9 有堤和無堤時波高及近底水平波動流速幅值沿程分布過程Fig.9 Spatial distribution of wave height and near-bottom horizontal fluctuation velocity amplitude in the presence and absence of a slope breakwater

為確定斜坡堤前護面塊體重量，《防波堤與護岸設計規范》（JTS 154—2018）［27］推薦按下式計算堤前最大波浪底流速：

式中：Vmax為斜坡堤前最大波浪底流速；H為設計波高；L為計算波長；d為堤前水深。在堤前半倍波長范圍內，由式（17）和無堤時區間最大波高計算的最大底流速為0.254 m/s，而數模計算的區間最大近底流速較之增大約70%；在礁坪上距斜坡堤腳0.87～0.76 倍波長、0.63～0.28 倍波長及接近堤腳處的流速大于規范值。因此在島礁地形上需注意加強堤前護底結構。

3 結 語

基于SPH 模型建立波浪數值水槽，對島礁地形上的波浪傳播變形進行了數值模擬，并與物理模型試驗數據進行對比，驗證了SPH 方法在模擬島礁地形上波浪傳播變形方面的可行性。比較分析了有、無斜坡堤情況下礁坪上的波浪破碎過程、流場形態及流速分布。

在無堤情況下，波浪在破碎前水平波動流速最大值自上而下逐漸減小，在礁坪上破碎后呈現出非線性特征，正向水平波動流速最大值自上而下逐漸減小，反向水平波動流速最大值自上而下逐漸增大。

有斜坡堤后，波浪破碎更加猛烈，并產生漩渦和紊動等復雜流體運動現象；波浪破碎點前移，同時波浪破碎后的紊動區前移，并且紊動更加劇烈；受反射波影響，礁坪上波高呈現為先增大、再減小、再增大，最后在防波堤堤腳處減小，近底流速變化趨勢則與之相反。在堤前半倍波長范圍內，有堤時反向近底流速峰值增幅較大，島礁上斜坡堤護底結構設計應關注這一現象。