不規(guī)則波作用下斜坡堤越浪的數(shù)值模擬

李東洋，張慶河，焦方騫

(天津大學(xué) 水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室，天津 300072)

斜坡堤是保護(hù)沿海地區(qū)免受潮、浪襲擊的重要水工建筑物，由于斜坡堤迎浪面較緩，波浪爬高往往較大，很容易發(fā)生越浪，因此準(zhǔn)確模擬斜坡堤越浪過程并確定越浪量，對于斜坡堤的設(shè)計具有重要意義。

早期對斜坡堤越浪的研究主要以物理模型試驗(yàn)為主，王聰?shù)萚1]在總結(jié)前人物理模型試驗(yàn)的基礎(chǔ)上，研究了不規(guī)則波作用下?lián)趵藟Ω叨鹊葘υ嚼肆康挠绊?。近年來，隨著計算機(jī)技術(shù)和并行計算的發(fā)展，利用數(shù)值模型研究斜坡堤越浪成為港口海岸工程學(xué)科關(guān)注的熱點(diǎn)。追蹤自由表面運(yùn)動的斜坡堤越浪數(shù)值模型主要分為兩大類，一類基于N-S方程，如Van Gant等[2]較早基于二維不可壓縮的N-S方程研究了波浪與寬肩臺可滲斜坡堤相互作用的水動力特性，Losada等[3]基于體積平均的RANS方程和VOF方法建立了COBRAS(Cornell Breaking Waves and Structure)模型，模擬分析了不規(guī)則波在可滲斜坡堤上的越浪過程。關(guān)大瑋[4]基于CFD軟件Flow-3D建立三維數(shù)值波浪水槽，模擬研究了規(guī)則波和不規(guī)則波作用下不可滲光滑斜坡堤的越浪過程。閆科諦等[5]基于OpenFOAM軟件建立數(shù)值波浪水槽，通過網(wǎng)格細(xì)化，刻畫出柵欄板凹槽的形狀，進(jìn)而模擬了規(guī)則波作用下非可滲堤心柵欄板護(hù)面的斜坡堤越浪過程。Xiang等[6]等通過將Fluidity和Y3D耦合的方法，模擬研究了規(guī)則波作用下扭王塊體護(hù)面斜坡堤周圍的水動力特性。另外一類模型為無網(wǎng)格方法，如葉曉文[7]基于不可壓縮SPH數(shù)學(xué)模型研究了規(guī)則波沿可滲斜坡堤的爬坡和越浪過程， Ren[8]基于弱可壓的光滑粒子流體力學(xué)方法(WCSPH)模擬了波浪與可滲潛堤的相互作用，并分析了堤體內(nèi)部的水動力特性。溫鴻杰[9]基于改進(jìn)的SPH方法，模擬研究了規(guī)則波在出水堤上的破碎和越浪過程。

總體上來看，目前關(guān)于斜坡堤越浪數(shù)值模擬的研究還主要集中于規(guī)則波，對護(hù)面塊體進(jìn)行三維全尺度模擬且考慮護(hù)面塊體墊層以下可滲特性的工作也比較少見，這使得利用數(shù)值模型解決斜坡堤越浪等實(shí)際工程問題仍受到一定的制約。因此，利用數(shù)值模型展現(xiàn)真實(shí)海況下斜坡堤越浪過程，需要進(jìn)一步探索。考慮到SPH類模型在壓力脈動等方面還存在有待解決的問題，本文將利用Jesus[10]和Lara等[11]基于OpenFoam開源軟件開發(fā)的ihFoam求解器，通過OpenFoam網(wǎng)格劃分工具snappyHexMesh精細(xì)刻畫網(wǎng)格的方法生成扭王字塊，模擬不規(guī)則波作用下扭王塊體護(hù)面斜坡堤越浪過程，以使斜坡堤越浪的數(shù)值模擬在解決工程問題時達(dá)到一定的實(shí)用性，為斜坡堤的設(shè)計提供依據(jù)。

1 控制方程與邊界條件

1.1 控制方程

本文通過將滲流運(yùn)動的Darcy-Forchheimer方程引入到兩相流求解器interFoam中，生成全新的考慮孔隙介質(zhì)流的ihFoam求解器，采用主動吸收式速度入口造波，通過求解體積平均的RANS方程，模擬水體和空氣兩相不可壓縮流體的運(yùn)動。紊流模型選用RNGk-ε模型，采用VOF法捕捉自由表面。模型控制方程如下

(1)

(2)

(3)

(4)

式中:D50為組成孔隙介質(zhì)的中值粒徑；KC=umT/(nD50)，其中um為最大震蕩速度；T為震蕩周期。α和β分別為Darcy-Forchheimer線性項(xiàng)系數(shù)和非線性項(xiàng)系數(shù)，2011年Lara等[13]提出了α和β關(guān)于孔隙率和中值粒徑的經(jīng)驗(yàn)公式如下

(5)

當(dāng)n=1時，即純水體情況時，式(2)轉(zhuǎn)化為標(biāo)準(zhǔn)的雷諾時均方程。

圖1 數(shù)值水槽和邊界條件示意圖Fig.1 Numerical wave flume and boundary condition diagram

1.2 邊界條件

數(shù)值水槽和邊界條件示意圖如圖1所示，左側(cè)入口設(shè)置為主動吸收造波邊界。右側(cè)出口設(shè)置為輻射消波邊界。數(shù)值水槽前后側(cè)邊界設(shè)置為可滑移固邊界，法向速度為0。底面及護(hù)面塊體設(shè)置為固壁邊界，數(shù)值邊界條件設(shè)置為無滑移邊界條件，壁面處的流速為0。

2 扭王塊體護(hù)面斜坡堤越浪模擬

為了驗(yàn)證數(shù)值模型的合理性，我們收集已有物理模型試驗(yàn)數(shù)據(jù)與數(shù)值模型結(jié)果進(jìn)行對比分析。

2.1 物理試驗(yàn)?zāi)Ｐ驮O(shè)置

圖2 物理模型試驗(yàn)尺寸示意圖[14]Fig.2 Physical model size diagram[14]

陳偉秋等[14]通過物模試驗(yàn)研究了不規(guī)則波作用下扭王塊體護(hù)面斜坡堤的平均越浪量。試驗(yàn)是在長175 m、寬1.2 m、高1.5 m的波浪水槽中進(jìn)行的，試驗(yàn)?zāi)Ｐ偷木唧w尺寸如圖2所示，斜坡堤左右兩側(cè)坡度均為1:1.5，堤頂高度為0.5 m，斜坡堤堤底寬為2.11 m，兩端各有坡度為1:2的拋石棱體作為護(hù)底，海側(cè)斜坡鋪設(shè)一層扭王字塊，單個扭王字塊的高度是6 cm。試驗(yàn)水深0.4 m，入射波為不規(guī)則波，采用JONSWAP譜進(jìn)行模擬，有效周期為2 s，有效波高分別為6.6 cm、7.1 cm、7.6 cm。

表1 多孔介質(zhì)物理參數(shù)表Tab.1 Physical parameter of porous media

模擬過程中除扭王塊體采用全尺度模擬外，斜坡堤堤體其他組成部分均作為孔隙介質(zhì)處理，各層物理參數(shù)見表1。

2.2 三維數(shù)值水槽不規(guī)則波驗(yàn)證

試驗(yàn)中采用JONSWAP譜進(jìn)行模擬，JONSWAP譜的表達(dá)形式為

(6)

(7)

(8)

(9)

式中：S(f)為波浪譜密度函數(shù);Hs為有效波高;fp為譜峰周期;f為波浪頻率;γ為譜峰因子，取3.3;σ是波浪密度的標(biāo)準(zhǔn)差，當(dāng)f≤fp時，σ=0.07，當(dāng)f>fp時，σ=0.09。

3-a 波面歷時圖 3-b 計算波譜與理論波譜對比圖圖3 Hs=7.1 cm情況下數(shù)值造波結(jié)果Fig.3 Hs=7.1 cm Simulation results of free surface

不規(guī)則波采用線性波疊加法生成[15]，根據(jù)等分頻率法將波浪能量分割成70份。因?yàn)閿?shù)值模擬的側(cè)壁采用的是可滑移邊界條件，避免了壁面對內(nèi)部流場的影響，所以可以利用寬度較窄的計算水槽反演試驗(yàn)結(jié)果，具體數(shù)值水槽設(shè)置如下：長15 m，高0.6 m，寬0.24 m，網(wǎng)格水平空間步長取L/100，為了準(zhǔn)確追蹤自由表面，在其靜水位上下2倍有效波高范圍內(nèi)取30層網(wǎng)格，上下網(wǎng)格漸變，寬度方向上取15層網(wǎng)格。數(shù)值模擬時間為220 s，取x=12 m(即堤腳)處測點(diǎn)進(jìn)行分析，以有效波高為7.1 cm為例，圖3-a顯示了其水位歷時曲線。取后100個波，即20～220 s的實(shí)測數(shù)據(jù)進(jìn)行快速傅里葉分析，得到不規(guī)則波浪的能量譜。圖3-b顯示了計算模擬波譜與理論波譜的對比情況，由圖可知兩者比較接近。模擬有效波高與理論有效波高誤差為2.4%，滿足《波浪模型試驗(yàn)規(guī)程》[16]要求。

圖4 AutoCAD三維建模中扭王字塊及扭王字塊排放示意圖 圖5 扭王塊體和胸墻經(jīng)snappyHexMesh處理后局部網(wǎng)格剖分示意圖Fig.4 AutoCAD three-dimensional model of accropode blocks Fig.5 Accropode blocks and parapet local mesh encryption by snappyHexMesh

2.3 扭王塊體全尺度模擬的網(wǎng)格刻畫

OpenFOAM提供了強(qiáng)大的網(wǎng)格生成工具，通過自帶的網(wǎng)格處理工具snappyHexMesh將扭王字塊和胸墻轉(zhuǎn)變成貼體網(wǎng)格模擬，可以逼近真實(shí)工況下斜坡堤上的波浪爬坡和越浪。

首先參照《防波堤設(shè)計與施工規(guī)范》[17]中規(guī)定的扭王字塊的尺寸，使用AutoCAD三維建模繪制出扭王字塊的外部輪廓。周雅等[18-19]通過物理模型試驗(yàn)研究得出扭王字塊護(hù)面塊體的排放型式對越浪量的影響，因此為了簡化模型，將扭王字塊規(guī)則排放在斜坡堤上，如圖4所示。圖5顯示了精細(xì)網(wǎng)格刻畫后的扭王塊體和胸墻局部網(wǎng)格形狀。由圖可以看出數(shù)模所采用的網(wǎng)格能較好描述扭王護(hù)面塊體和胸墻的外部形狀。

6-a越浪量歷時圖6-b累計越浪量圖圖7 平均越浪量計算值與試驗(yàn)值比較圖6 Hs=7．1cm情況下越浪量計算結(jié)果Fig．6Hs=7．1cmSimulationresultsofovertoppingdischargeFig．7Comparisonofnumericalaveragedischargeandexperimentalaveragedischarge

2.4 越浪模擬與試驗(yàn)結(jié)果的比較

8-a t=174.36 s 9-a t=174.36 s

圖6顯示了有效波高為7.1 cm情況下的波浪越浪量歷時圖6-a和累計越浪量圖6-b，在有效周期一定的情況下，越浪量隨有效波高的增大而明顯增大。取20～220 s(100個波)時間內(nèi)的越浪量計算平均值以得到平均越浪量。圖7顯示了平均越浪量數(shù)值模擬結(jié)果與試驗(yàn)值的比較，由圖可知斜坡堤越浪的數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果吻合較好，表明利用snappyHexMesh網(wǎng)格處理工具生成貼體網(wǎng)格的方法對扭王護(hù)面塊體進(jìn)行全尺度模擬，對墊層以下塊石利用滲流模型模擬斜坡堤越浪可以獲得較好的效果。

2.5 越浪過程和堤前附近流場分析

8-b t=174.56 s 9-b t=174.56 s

8-c t=174.76 s 9-c t=174.76 s圖8 最大單波越浪過程三維示意圖 圖9 堤前附近流場圖 Fig.8 3D process of the largest Fig.9 Variation diagram of single-wave overtopping velocity fields

圖8顯示了有效波高為7.6 cm，在174.5 s前后發(fā)生越浪時波面的三維過程，圖9顯示了Y=0.1 m，XZ剖面扭王字塊附近的流場變化。從圖中可以看出，越浪發(fā)生前，波浪沿著坡面向上爬升，在向上水體的帶動下，水面抬高，波浪拍擊胸墻表面，直至越過堤頂形成越浪。在整個越浪過程中，由于扭王字塊護(hù)面的作用，堤內(nèi)的滲流流速相對很小。

從越浪過程、越浪波面和流速變化過程看，本文所建立的模型較準(zhǔn)確地反映了波浪與扭王塊護(hù)面斜坡堤的相互作用，具備了為斜坡堤設(shè)計優(yōu)化提供支持的能力。

表2模型與原型波要素

Tab.2 Wave elements under experimental model and prototype size

類別有效波高(m)有效周期(s)水深(m)實(shí)驗(yàn)室尺度0．07120．40原型尺度0．071×25=1．7752×251/2=100．40×25=10．00

2.6 原型工況的數(shù)值模擬

現(xiàn)階段對于原型尺寸下的波浪與斜坡堤相互作用的研究，大多是通過比尺縮放，進(jìn)行實(shí)驗(yàn)室物模試驗(yàn)，再將試驗(yàn)值根據(jù)《波浪模型試驗(yàn)規(guī)程》換算成原型，而直接研究原型尺寸下斜坡堤越浪的工作還比較少。因此，為進(jìn)一步論證模型可否應(yīng)用于描述現(xiàn)場尺度條件下的斜坡堤越浪，首先依據(jù)《波浪模型試驗(yàn)規(guī)程》[16]將文獻(xiàn)[14]中的物理模型試驗(yàn)按照1∶25的比尺放大成原型工況，通過模擬該尺度下的斜坡堤越浪，并與文獻(xiàn)[14]按照重力相似計算得出的越浪量進(jìn)行比較，從而驗(yàn)證模型是否可用于模擬計算現(xiàn)場尺度下的斜坡堤越浪。以模型有效周期為2 s，有效波高0.071 m為例，按照重力相似準(zhǔn)則確定原型波要素如表2所示。

圖10 原型工況下累計越浪量Fig.10 The accumulated overtopping discharge under prototype size

原型模擬數(shù)值水槽長×寬×高=375 m×6 m×15 m，依然按照上述方法建模。圖10顯示了原型尺度下累計越浪量的數(shù)值模擬結(jié)果，單個越浪發(fā)生時其實(shí)歷時很短，可以看成是一個瞬時過程，圖中折線每次的抬升就代表發(fā)生了一次越浪，“臺階”的高度即為單個越浪發(fā)生時的越浪量，從而可以得出，單波最大越浪量發(fā)生在950 s左右，單波最大越浪量約為1.3 m3/m。取100～1 100 s(100個波)時間內(nèi)的越浪量計算平均值得到平均越浪量，原型尺度下數(shù)值計算的越浪量為7.89×10-3m3/m/s。根據(jù)模型試驗(yàn)越浪量折算到原型的平均越浪量為6.68×10-5×251.5=8.36×10-3m3/m/s，與原型數(shù)值模擬結(jié)果較為接近。合理的模型試驗(yàn)可以反映原型結(jié)果是比較公認(rèn)的，因此，從原型直接模擬結(jié)果與模型試驗(yàn)折算到原型的越浪量比較吻合來看，利用OpenFOAM直接模擬原型工況下的斜坡堤越浪是可行的，并推薦采用原型尺寸進(jìn)行數(shù)值研究，進(jìn)而優(yōu)化確定設(shè)計階段中斜坡堤的斷面型式和頂高程等。

3 結(jié)語

本文基于OpenFOAM建立數(shù)值波浪水槽，控制方程采用了體積平均的RANS方程和包含非線性項(xiàng)的滲流方程描述斜坡堤內(nèi)部的水體滲流流動。利用AutoCAD三維建模繪制出扭王塊體的外部輪廓，利用snappyHexMesh網(wǎng)格劃分工具對斜坡堤扭王護(hù)面塊體進(jìn)行全尺度網(wǎng)格精細(xì)劃分。采用建立的模型模擬研究了正向入射不規(guī)則波與扭王塊體護(hù)面斜坡堤的相互作用，模擬結(jié)果表明，平均越浪量的計算值與實(shí)驗(yàn)結(jié)果吻合較好。通過將實(shí)驗(yàn)室尺度放大成原型尺度，模擬研究了原型條件下的斜坡堤越浪，從計算結(jié)果可以看出，平均越浪量的計算值與模型通過重力相似準(zhǔn)則折算成原型的結(jié)果較為吻合。從實(shí)驗(yàn)室尺度和原型尺度模擬結(jié)果看，對護(hù)面塊體進(jìn)行全尺度模擬的數(shù)值波浪水槽目前已可以較為合理地描述復(fù)雜護(hù)面塊體斜坡堤的越浪過程，可以應(yīng)用于設(shè)計階段對斜坡堤斷面型式和頂高程的優(yōu)化確定。

需要指出的是，本文僅給出了波浪正向入射情況下斜坡堤的越浪模擬結(jié)果，有關(guān)護(hù)面塊體本身的穩(wěn)定性以及斜向入射波與斜坡堤的相互作用有待進(jìn)一步研究。

[1]王聰, 陳國平, 嚴(yán)士常, 等. 不規(guī)則波作用下斜坡堤越浪量試驗(yàn)研究[J]. 水運(yùn)工程, 2017(2): 49-52.

WANG C, CHEN G P, YAN S C,et al. Experimental study on wave overtopping of sloping seawall under irregular wave action[J]. Port & Waterway Engineering, 2017(2): 49-52.

[2]Van Gent M R A. Wave action on and in permeable structures[D]. The Netherlands: Delft University, 1995.

[3]Losada I J, Lara J L, Guanche R, et al. Numerical analysis of wave overtopping of rubble mound breakwaters[J]. Coastal Engineering, 2008, 55(1): 47-62.

[4]關(guān)大瑋. 規(guī)則波與不規(guī)則波的海堤越浪數(shù)值模擬[D]. 廣州：華南理工大學(xué), 2016.

[5]閆科諦, 張慶河. 柵欄板護(hù)面斜坡堤越浪數(shù)值模擬研究[J]. 中國港灣建設(shè), 2016, 36(2):16-19.

YAN K D, ZHANG Q H. Numerical simulation of overtopping on sloping dike with fence panels[J]. China Harbour Engineering, 2016, 36(2):16-19.

[6]Xiang J, Latham J P, Vire A, et al. Coupled fluidity/Y3D technology and simulation tools for numerical breakwater modelling[J]. Coastal Engineering Proceedings, 2012, 1(33):66.

[7]葉曉文. 斜坡式防波堤越浪過程的SPH模擬[D]. 大連：大連理工大學(xué), 2011.

[8]Ren B, Wen H, Dong P, et al. Improved SPH simulation of wave motions and turbulent flows through porous media[J]. Coastal Engineering, 2016, 107: 14-27.

[9]溫鴻杰. 非線性波浪與可滲結(jié)構(gòu)相互作用的SPH模型[D]. 大連：大連理工大學(xué), 2016.

[10]Jesus M D, Lara J L, Losada I J. Three-dimensional interaction of waves and porous coastal structures. Part I: Numerical model formulation[J]. Coastal Engineering, 2012, 64:57-72.

[11]Lara J L, del Jesus M, Losada I J. Three-dimensional interaction of waves and porous coastal structures. Part II: Experimental validation[J]. Coastal Engineering, 2012, 64: 26-46.

[12]Higuera P, Lara J L, Losada I J. Three-dimensional interaction of waves and porous coastal structures using OpenFOAM?. Part I: Formulation and validation[J]. Coastal Engineering, 2014, 83: 243-258.

[13]Lara J L, Losada I J, Maza M, et al. Breaking solitary wave evolution over a porous underwater step[J]. Coastal Engineering, 2011, 58(9):837-850.

[14]陳偉秋, 陳凌彥, 王登婷,等. 斜坡堤后坡砌石護(hù)面穩(wěn)定厚度的模型試驗(yàn)研究[J]. 水運(yùn)工程, 2016(6): 93-98.

CHEN W Q, CHEN L Y, WANG D T, et al. Experimental study on thickness of inner block revetment of sloped seawall[J]. Port & Waterway Engineering, 2016(6): 93-98.

[15]俞聿修, 柳淑學(xué). 隨機(jī)波浪及其工程應(yīng)用[M]. 大連：大連理工大學(xué)出版社, 2011.

[16]JTJ/T234-2001，波浪模型試驗(yàn)規(guī)程[S].

[17]JTS154-1-2011, 防波堤設(shè)計與施工規(guī)范[S].

[18]周雅, 林登榮, 李慶銀, 等. 不規(guī)則波作用下斜坡堤越浪量試驗(yàn)研究[J]. 水道港口, 2016(4)：331-335.

ZHOU Y, LIN D R, LI Q Y, et al. Experimental research of wave overtopping on sloping dike under irregular waves[J]. Journal of Waterway and Harbor, 2016(4): 331-335.

[19]鐘雄華，陳國平，嚴(yán)士常，等.不同擺放方式扭王字塊體穩(wěn)定性研究[J].水道港口，2016(5)：479-483.

ZHONG X H, CHEN G P, YAN S C,et al. A stability analysis of different stack types of accropode[J]. Journal of waterway and harbor, 2016(5):479-483.