孤立波低頂海堤越浪的數值模擬研究

劉竹琴，殷銘簡，趙西增，2，羅 敏

（1.浙江大學 海洋學院，浙江 舟山 316021；2.浙江海洋大學 海洋工程裝備學院，浙江 舟山 316022）

隨著全球氣候變化，災害性海浪（臺風浪、海嘯等）發生的頻率和強度呈上升趨勢。圖1為強臺風“煙花”登錄浙江溫嶺產生的災害性臺風浪［1］，其引發的越浪對海岸附近人員安全造成嚴重威脅，誘發的海水倒灌對海岸及岸后結構物造成嚴重破壞，使得海防工程承受越來越大的風險。

圖1 強臺風“煙花”登陸浙江溫嶺[1]（圖片來源：https://new.qq.com/rain/a/20210723A02VJ700）Fig.1 Severe Typhoon In-Fa landing in Wenling, Zhejiang[1]（Source：https://new.qq.com/rain/a/20210723A02VJ700）

浙江沿海地區普遍采用海堤來防止海浪的破壞性作用（如越浪和劇烈波浪撞擊）。在常規海浪情況下，海堤可以對海岸起到很好的防護作用。然而在臺風期間，風暴潮導致水面上升，使得現役海堤堤頂高度減小，轉變為低頂結構（即按EurOtop（2018）［2］定義，相對超高為0 ≤Rc/H0< 1.5的狀態），無法有效防止越浪；此外，劇烈海浪撞擊海堤，產生巨大荷載，對海堤結構、附近基礎設施和人員安全造成嚴重威脅。

孤立波被廣泛應用于模擬長波等災害性海浪在淺水中的運動［3］。孤立波經過海堤產生的越浪量和波浪力是目前評判越浪作用機制的主要參數。對于孤立波越浪量經驗公式的研究，Ozhan 和Yalgmer［4］試驗研究孤立波經過坡面較陡的三角海堤的爬高和越浪量，并基于堰流比擬的方法提出了越浪量經驗公式，但不具有普適性；目前國內外普遍認可的經驗公式為Baldock等［5］通過大量試驗研究總結得到的孤立波爬高和越浪量經驗公式，研究表明越浪量與爬高具有較高相關性，因此在計算越浪量時需采用爬高數據；隨后，張金牛［6］試驗研究孤立波經過梯形海堤越浪量，提出的經驗公式依賴于相對波高（0.144 ≤H0/h0≤ 0.607）和相對超高（0.165 ≤Rc/H0≤ 2.649）；此外，Lee 等［7］通過物理模型試驗得到孤立波經過垂直護岸和斜坡護岸的越浪量經驗公式；魏斐斐等［8］采用FLUENT 軟件，模擬孤立波經過斜坡堤的越浪，并總結出越浪量沿堤頂寬度衰減的經驗公式。越浪量是海堤工程設計中的一個關鍵因素，國內外學者對該物理量開展了大量研究，并提出如下關于孤立波越浪量的經驗公式：

1）Baldock 等［5］通過物理模型試驗，綜合分析試驗數據以及前人研究成果，認為孤立波的越浪量受爬高和堤頂超高的影響，得到計算孤立波爬高和越浪量的經驗公式：

式中：q為孤立波的單寬越浪量；H為波高；d為水深；R為孤立波的爬高；Rc為堤頂超高；q0表示孤立波在傳播時相對于靜水面的孤立波水體體積，其計算公式見式（17）。該公式在計算越浪量時需采用爬高數據。

2）張金牛［6］通過物理模型試驗，研究不同水深和波高下孤立波經過梯形海堤的越浪量，認為越浪量與相對波高和相對超高有關，采用Baldock 等［5］方法對其進行無量綱化，得到經驗公式，該公式適用范圍為0.165 ≤Rc/H≤ 2.649。

Baldock等［5］提出的經驗公式需先計算孤立波的爬高，再通過爬高計算越浪量，較為繁瑣，不能通過水深和波高數據進行直接求解；張金牛［6］提出的公式相對超高適用范圍為0.165 ≤Rc/H≤ 2.649，不適用于0 ≤Rc/H≤ 0.164。且我國《海堤工程設計規范》（GB/T 51015—2014）［9］以及EurOtop（2018）規范［2］缺少計算孤立波越浪量經驗公式。因此，文中研究旨在填補孤立波在低頂海堤下越浪量經驗公式的空白。

對于波浪力的研究，Hsiao 和Lin［10］研究3 種相對波高的孤立波經過位于斜坡海岸上海堤的越浪形態及其在海堤兩側產生的波浪力，并分析孤立波在海堤后方的渦流運動與卷入氣體的關系，試驗數據可用于校準數值模型。其中，Tsung 等［11］以Hsiao 和Lin［10］的研究為基準算例，采用高階Boussinesq 方程研究不同水深下孤立波經過不透水海堤產生的波浪力和越浪量，并討論斜坡坡度以及波浪非線性作用的影響；此外，Luo等［12］也以Hsiao 和Lin［10］的研究為基準算例，采用CPM 模型模擬孤立波經過海堤時產生的波浪力和越浪量，并討論海堤前坡坡度對越浪的影響。但未具體討論孤立波作用在海堤兩側波浪力與相關參數的變化規律，對此需開展進一步研究。

基于上述背景，文中采用已被證實能夠有效模擬波浪與結構物作用的開源程序OpenFOAM 中interFoam求解器［13］，以浙江沿海普遍采用的斜坡堤海堤斷面型式為原型，考慮一種代表性長波即孤立波，研究孤立波在不同波高、水深和堤頂超高下的越浪過程，揭示越浪機理，總結出適用于低頂海堤的孤立波越浪量經驗公式。并探究防浪墻高度對越浪特征的影響以及防浪墻所受波浪力的變化，綜合考慮防浪墻減少越浪以及自身所受波浪力，給出建議的防浪墻高度。

1 數值模型理論和方法

采用OpenFOAM 中的兩相流求解器interFoam，求解雷諾時均連續性方程和Navier-Stokes 方程，采用有限體積法（finite volume method）進行數值離散，模擬水體和空氣兩相不可壓縮流體的運動。

式中：μ為動力黏度，k為湍流動能，δij為克羅內克符號，υt為湍流運動黏度。

采用流體體積（VOF）法追蹤自由液面。定義體積分數α，α= 1 表示該單元充滿水，α= 0 表示該單元充滿空氣，α介于0和1之間表示單元中既有空氣又有水。通過求解α的輸運方程，即式（8）來追蹤自由液面：

為修正尖銳表面，且限制α的取值在0和1之間，引入人工壓縮項對式（8）進行修正：

式中：uc為壓縮速度，uc，i= min(cα|ui|，max|ui|)，其中cα是控制交界面壓縮程度的參數，這里取值為1［13］。該人工壓縮項是守恒的，且在全是水或空氣時取值為0，即只在水體與空氣的交界面處生效。

基于體積分數α的加權值計算得到每個單元內流體的密度ρ：

式中：ρwater和ρair分別表示水和空氣的密度。

文中模擬采用Devolder等［14］建立的浮力修正SST k?ω湍流模型，基本方程為：

式中：k為湍動能，ω為耗散率，υ為流體的運動黏度，S為流體的平均應變率；β*= 0.09；a1= 0.31；F1和F2為計算σω、σk、β和γ所用到的參數，取值詳見Devolder等［14］研究。

2 數值模型驗證

2.1 模型設置

通過模擬具有試驗［10］和數值基準數據的孤立波與海堤相互作用的算例來驗證所采用的數值模型的準確性。二維數值波浪水槽的計算域如圖2所示，水槽長22 m，高0.4 m。水槽左側為造波區，在距離造波邊界7 m處有一坡度比為1∶20的斜坡，在距離造波邊界10.6 m處的斜坡上有一個海堤。在海堤前、后坡不同位置處放置壓強測點用于記錄海堤上的波浪動壓強，并在數值水槽代表位置放置浪高儀用于測量波高。用于測量波高的浪高儀和動壓強的測點的位置如圖2所示，其詳細位置如表1和表2所示。選取Hsiao和Lin［10］研究中水深h0=0.256 m、波高H0=0.058 9 m（波陡ε=H0/h0=0.23）的孤立波工況進行模擬。

表1 浪高儀位置Tab.1 Positions of wave gages

表2 壓強測點坐標位置Tab.2 Positions of pressure sensors

圖2 計算域尺寸、浪高儀和壓強測點位置以及邊界條件Fig.2 Dimensions of computational domain, positions of wave gauges and pressure sensors and boundary conditions

關于計算邊界，計算域底部（包括結構物表面）和出口采用無滑移邊界，計算域頂部采用零梯度邊界，計算域入口采用速度邊界法產生孤立波。文中采用一階孤立波公式，波面方程η為：

2.2 收斂性分析

采用3套網格進行網格收斂性驗證，分別為細網格、中網格、粗網格，3套網格參數信息見表3，均采用自適應時間步，庫朗數上限為0.25。不同網格尺寸下的g22處波高歷時曲線如圖3所示，p7處動壓強歷時曲線如圖4所示。模擬結果表明，粗網格下的波高和壓強結果明顯比其他兩組數值結果偏小，中網格和細網格下的數值結果趨于一致。對于中網格，進行時間步長收斂性分析，圖5 和圖6 分別為庫朗數上限為0.10、0.25、0.50 下g22 處波高歷時曲線和p7 處動壓強歷時曲線結果對比，由圖可知庫朗數上限0.10 與0.25 結果趨于一致，庫朗數上限0.50 模擬結果偏小。綜合考慮數值模擬精度和效率，選取中網格尺寸以及庫朗數上限0.25開展模擬研究。

表3 收斂性驗證網格信息Tab.3 Convergence verification grid information

圖3 不同網格尺寸下g22處波高歷時曲線Fig.3 Wave elevations at g22 under different grid sizes

圖4 不同網格尺寸下p7動壓強歷時曲線Fig.4 Dynamic pressures at p7 of the seawall under different grid sizes

圖5 不同庫朗數下g22處波高歷時曲線Fig.5 Wave elevations at g22 under different Courant numbers

圖6 不同庫朗數下p7動壓強歷時曲線Fig.6 Dynamic pressures at p7 of the seawall under different Courant numbers

2.3 數值模擬與試驗結果對比分析

圖7 所示為波高歷時曲線數值模擬結果與試驗結果對比，波高峰值的平均相對誤差為3%；圖8 所示為孤立波對海堤的動壓強數值模擬結果與試驗結果對比，動壓強峰值的平均相對誤差為4.9%；圖9 所示為孤立波對海堤兩側的動壓力數值模擬結果與試驗結果對比，動壓力峰值的平均相對誤差為5.1%。總的來說，OpenFOAM 模擬結果與試驗結果吻合良好。但是對于后坡沖擊壓強和壓力，OpenFOAM 峰值略小于試驗峰值。這是由于數值模擬基于二維不可壓縮模型，難以準確模擬該過程中的水氣摻混現象，導致計算的壓強和壓力與試驗結果具有偏差?？傮w來說，OpenFOAM能夠較為準確模擬孤立波越浪。

圖7 H0=0.256 m, h0=0.058 9 m下典型位置的波高歷時曲線Fig.7 Wave elevations at typical positions of H0=0.256 m, h0=0.058 9 m

圖9 H0=0.256 m, h0=0.058 9 m下海堤兩側動壓力Fig.9 Dynamic forces on the seawall of H0=0.256 m, h0=0.058 9 m

3 浙江省某海堤孤立波越浪的數值模型設置

以我國浙江省沿海普遍采用的斜坡堤海堤斷面型式為原型，對浙江省海岸的越浪特性開展研究。綜合考慮計算效率和Froude相似準則，采用縮尺比為λ=16的海堤進行數值模擬。圖10為數值波浪水槽示意，在距離造波邊界10 m處有一坡比為1∶20的斜坡，斜坡長度為5 m；斜坡后方連接海堤，海堤前坡坡比為1∶3，長度為0.6 m；海堤堤頂高度為0.45 m，寬度為0.325 m；海堤后坡坡比為1∶2，長度為0.4 m；海堤后方為平底地形，高度為0.25 m。

圖10 中國浙江省典型海堤數值波浪水槽示意Fig.10 Schematic diagram of the numerical wave flume with the typical seawall of Zhejiang Province, China

浙江省近岸區域的水深為6.4～7.2 m、有效波高為0.96～1.92 m。根據縮尺比，研究0.40、0.42、0.44 和0.45 m 四組水深以及0.06～0.12 m 七組波高，共28組工況，如表4所示，其中qmax為最大越浪率，q*為無量綱單寬越浪量。

表4 28組工況下的孤立波參數及最大越浪率和無量綱單寬越浪量Tab.4 Solitary wave parameters, maximum overtopping rate and dimensionless unit width overtopping volume of 28 cases studied

海堤的堤頂超高Rc指堤頂到靜止水面的距離，可通過改變水深實現不同的堤頂超高，堤頂超高與波高的比值為海堤相對超高Rc/H0。

4 結果與討論

4.1 越浪形態

圖11所示為h0=0.40 m（Rc=0.05 m）、H0=0.12 m工況，孤立波經過海堤的代表性時刻波面形態。孤立波峰在8.84 s時達到海堤堤趾，由于海堤前坡水深急劇變淺，波峰前側波高增大（圖11（a））；波浪沿著前坡爬升至海堤堤頂，未發生破碎（圖11（b）），并直接越過海堤堤頂（圖11（c）），而后沿著海堤后坡流下（圖11（d）），并沖擊海堤后坡（圖11（e））以及海堤后方的建筑物和設施（圖11（f））。

圖11 工況7中代表性時刻的孤立波爬升和越浪形態Fig.11 Profiles of solitary wave run-up and overtopping at typical time instants in Case 7

圖12 所示為h0=0.45 m（與堤頂高程相同）、H0=0.12 m 工況，孤立波經過海堤的代表性時刻波面形態變化。孤立波峰在8.68 s時到達海堤堤趾，由于水深變淺，波高增大，在海堤前坡未發生破碎（圖12（a））。越浪水體到達海堤堤頂時，直接越過海堤堤頂（圖12（b））；隨著越浪水體的進一步傳播，波浪越過海堤頂部沿著后坡流下（圖12（c），圖12（d）），直接沖擊海堤后坡（圖12（e）），對海堤后的建筑物和設施造成嚴重破壞，并引發內陸洪水（圖12（f））。相較于h0=0.40 m，當水深與堤頂高程相同時，波浪到達時間提前，越過堤頂所需時間減少，越過的水體增加，對海堤后方造成的危害加大。

圖12 工況28中代表性時刻的孤立波爬升和越浪形態Fig.12 Profiles of solitary wave run-up and overtopping at typical time instants in Case 28

4.2 越浪量

在數值模擬中，統計孤立波經過海堤堤頂與后坡交界處（x=15.925 m）在不同時間下的單寬越浪率（m3·s?1·m?1），對時間積分得到單寬越浪量q（m3/m）。采用Baldock 等［5］的公式形式，引入孤立波在靜水面以上的單寬水體體積q0，最終計算得出單寬無量綱越浪量q*，對式（16）積分可得q0為：

式中：h0為水深，H0為孤立波波高。

單寬無量綱越浪量q*的計算式為：

采用式（18），計算得到文中所有工況下的無量綱單寬越浪量q*。表4所示為各工況下的最大越浪量qmax和無量綱單寬越浪量q*。

圖13 所示為最大越浪量與水深、波高的關系。由圖13 可知，隨著水深和波高的增大，最大越浪量顯著增大。當水深與海堤堤頂高度接近或齊平時，即使很小的波也會產生較大越浪量，對海堤上及堤后的人員和建筑物造成嚴重威脅。

圖13 不同水深和波高下的最大越浪量Fig.13 The maximum overtopping volume at various water depths and wave heights

前人關于越浪量經驗公式的研究認為，影響越浪量的因素主要分為兩類：第一類為波浪要素，包括水深、入射波高、波周期、波陡等；第二類為海堤幾何尺寸，包括海堤堤頂高度、堤頂寬度等。

文中模擬的海堤表面光滑且幾何尺寸不變，對越浪量進行無量綱化處理，得到越浪量與各影響因素之間的函數關系式為：

根據π 定理，由量綱分析得到如下的無因次函數關系：

通過對前人的研究成果進行對比分析，同時考慮到文中模擬為堤頂超高很小或為0情況，提出的孤立波越浪計算公式采用張金牛［6］提出的結構形式為：

式中：a，b，c均為經驗參數。采用多元線性回歸方法對數值模擬所得的越浪量數據進行擬合，得到如下預測公式：

圖14為孤立波作用下海堤無量綱越浪量的數值模擬結果與公式預測結果的比較，由圖可知二者結果吻合較好，擬合的均方根誤差僅為0.015 3，表明公式預測結果與數值模擬結果具有較高的一致性。Baldock等［5］提出采用爬高計算越浪量，在計算h=0.40 m、0.42 m 時越浪量為負值，與實際不符，式（22）的計算公式更接近低頂情況下的越浪量。

圖14 無量綱單寬越浪量模擬結果與公式預測結果對比Fig.14 Comparison between the simulated dimensionless unit width overtopping volume and the predicted results by the proposed formula

4.3 越浪引起的波浪力

孤立波拍打在海堤上時，產生相當大的波浪力，對海堤的結構安全產生威脅。圖15為水深0.40 m、不同波高下，海堤前坡和后坡所受的單寬水平和垂向波浪力。結果表明，對于前坡，其所受的垂向力為水平力的3 倍，這是由前坡的坡度決定。隨著波高的增加，前坡受到的波浪力顯著增大，其中水平力從波高0.06 m 時的0.165 kN/m 增加到波高0.12 m 時的0.243 kN/m，垂向力從波高0.06 m 時的0.495 kN/m 增加到波高0.12 m時的0.728 kN/m，增長了47%。由于海堤前坡的坡度較緩，前坡所受水平力較小，具有較好的抗傾覆能力。此外，這種緩坡海堤對孤立波的阻力也較小，引起的波浪反射較少。因此，波浪會以較大的速度在前坡爬升并越過堤頂，對海堤后方的建筑物和設施產生較大的波浪力，造成更大的破壞。對于后坡，由于孤立波經過海堤時，其所受的波浪力比前坡小，后坡所受的垂向力為水平力的2 倍。隨著波高的增加，后坡受到的波浪力也增大，增長幅度為177%，遠大于前坡，后坡結構更易受到波浪力的影響。

圖15 水深0.40 m，不同波高下海堤兩側波浪力Fig.15 At water depth of 0.40 m, wave forces on both sides of the seawall under different wave heights

圖16 為不同水深和波高下的最大波浪力。由圖16 可知，在相同波高下，隨著水深的增加，海堤前坡所受的波浪力無顯著差異，海堤后坡所受的波浪力隨水深變化明顯。這是由于海堤后坡所受波浪力由越過堤頂的波浪產生，水深的增大導致越過堤頂的波浪增加，對后坡產生的波浪力也越大。

圖16 不同水深、波高下的最大波浪力Fig.16 The maximum wave forces on the seawall at different water depths and wave heights

4.4 防浪墻高度對越浪的影響

由于海平面上升，浙江沿海海堤堤頂高度不滿足海堤建設規范要求的可容忍限度越浪量。為了應對海平面上升以及極端海況下的越浪問題，一種方式是在海堤頂部建造防浪墻（通常是在局部海岸建造），以增加堤頂超高，在增加少量成本的情況下提高安全性。因此，文中提出在海堤堤頂向海一側（x= 15.6 m）設置一個防浪墻（如圖17所示），研究防浪墻寬度dw=0.01 m、高度hw分別為0.10、0.12、0.14、0.16、0.18、0.20 m在最不利工況28下的越浪特性。

圖17 帶防浪墻的海堤斷面示意Fig.17 Schematic diagram of the seawall with a crown wall

圖18為孤立波經過防浪墻高度為0.20 m 的波面形態演變。當孤立波爬升至海堤堤頂時，由于防浪墻的存在，無法直接越過堤頂，而是沖擊墻體并沿著墻面爬升（圖18（a））；當波浪達到防浪墻頂部時（圖18（b）），水體繼續往上升，超過防浪墻頂部（圖18（c））；當達到極限高度時，水體坍塌，少量水體回落至前坡，大部分水體從防浪墻后側落下（圖18（e））；越過防浪墻的水體在重力作用下產生堰流沖擊海堤堤頂（圖18（f））；堰流沖向堤頂破碎，分成兩股流（圖18（g）），一股沿著后坡流下，另一股流向防浪墻后側墻面（圖18（h））。

圖18 防浪墻高度0.20 m的波面形態演變Fig.18 Wave profiles with crown wall height of 0.20 m

波浪越過墻體的越浪量及其對防浪墻的波浪力與防浪墻高度有關。使用國內外學者普遍認可的方式［15-17］對防浪墻高度以及最大波浪力進行無量綱化，即防浪墻高度Hw除以波高H0得到無量綱防浪墻高度，單寬最大波浪力除以ρgH0dwhw得到無量綱化單寬最大波浪力F*。圖19為不同無量綱防浪墻高度下的單寬越浪量和最大單寬越浪波浪力，其中防浪墻高度Hw*=0 表示海堤上方無防浪墻的情況。相較于無防浪墻的海堤，防浪墻的存在能夠有效減少孤立波的越浪，極大減少越浪對海堤后方建筑物及人員的破壞作用。具體來說，當防浪墻高度為0.10～0.14 m 時越浪量線性減小，進一步增加防浪墻高度越浪量減小趨勢減緩。然而，防浪墻所受的波浪力隨其高度呈現先減小后增大而后又減小的趨勢：Hw*為0.83時受到最大波浪力q*=0.133，Hw*=1.00 受到最小波浪力q*=0.115；進一步增加Hw*受到的波浪力增大，在Hw*=1.33 達到最大值q*=0.124；而后波浪力隨Hw*減小。綜合考慮防浪墻減少越浪效果以及自身所受波浪力，針對文中采用的海堤截面和波浪條件，建議為1.00。

5 結 語

采用OpenFOAM 建立數值水槽，模擬研究孤立波的低頂海堤越浪特征以及越浪在海堤兩側產生的波浪力，同時探究海堤頂部防浪墻高度對越浪的影響。

所有工況下的波浪均未在海堤前坡發生破碎，孤立波越過堤頂，沿海堤后坡流下，直接沖擊堤后建筑物和設施。堤頂超高減小導致更為劇烈的越浪，這意味著在全球海平面上升的背景下，沿海城市面臨著更大的海堤越浪和海水倒灌風險。針對低頂海堤孤立波越浪量經驗公式的不足，文中提出適用于堤頂超高小或為0的孤立波越浪量經驗公式，公式預測結果較為準確。

在海堤頂部設置防浪墻增加海堤堤頂超高，可有效減少越浪，但防浪墻所受的波浪力也增大。綜合考慮防浪墻減少越浪效果以及自身所受波浪力，針對文中采用的海堤截面和波浪條件，建議無量綱防浪墻高度為1.00。