不規則波在孔隙介質礁坪上傳播過程的波高和增水變化數值研究

何棟彬，馬玉祥*，董國海

( 1. 大連理工大學 海岸和近海工程國家重點實驗室，遼寧 大連 116024)

1 引言

自然條件的島礁地形上孔隙介質相當常見，塊狀的巖石、珊瑚骨骼結構乃至植被，都可以視為孔隙介質；在人工島礁或者防波堤等工程結構中，也往往采用孔隙介質來減弱波浪的沖刷，從而起到保護海工建筑物的作用。然而孔隙介質的影響在以往的數值模型中往往被忽略，或者使用一個包含拖曳力系數Cd和速度平方項的關系式來進行處理[1-3]，這種方式較為簡單但比較粗糙。近些年來，使用體積平均的雷諾時間平均（Volume-Averaged Reynold-Averaged Navier-Stokes, VARANS）方程來處理孔隙介質中的流動問題[4]，由于其沒有引入靜壓假定、緩坡假定等近似假設，并完整保留了垂向動量方程，對各類地形水深條件不存在限制，適用于復雜水波運動，正受到越來越多的關注。

孔隙介質中存在流體時，被認為是一種包含固體和液體的兩相介質，通常假設為剛性的，內部空腔互相連通[5]?？紫督橘|的一大特征是其內部空腔的形狀、大小不規則，且分布無規律，所以模擬孔隙內部每一個流體質點的運動情況是不現實的，也是沒有必要的。VARANS模型的解決思路是對包含孔隙介質的流動區域在空間上進行體積平均，通過平均化的光滑效果，將孔隙介質的空間不均勻性所導致的小尺度波動過濾掉，類似于雷諾平均方法處理湍流問題。VARANS方程與雷諾時間平均（Reynold-Averaged Navier-Stokes, RANS）方程通過孔隙率（孔隙介質內部空腔的體積占比）相聯系，可采用同一形式表達，統一了孔隙介質內部和純流體區域的流動方程形式，即在孔隙介質內部，流動方程為VARANS方程，但在孔隙介質外部，VARANS方程自動轉換為RANS方程。經過體積平均后的RANS方程需要有相關的封閉項對方程進行封閉，經過多年的發展，現在基本采用經過擴展的Forchheimer關系式，該表達式首先由Forchheimer提出，即在原始的Darcy公式的基礎上，加入考慮湍流運動的阻力項[6]，后續Polubarinova-Kochina[7]又加入了包括附加質量系數的時間偏導數項來考慮流體的非定常流動效應，其中的各項系數經過大量的研究和校正[8]，目前取值已趨于固定，并得到了較為廣泛的驗證[9-12]。另外，對于VARANS方程中的湍流耗散項，體積平均后也需要對相應的湍流模型進行處理，目前由Nakayama和Kuwahara[13]建立的體積平均的k-ε模型是最常用的。

在多種水波數值模型中，非靜壓模型由于其將自由表面處理為單值函數，簡化了追蹤自由表面的難度并提升了計算效率，同時保持了很高的精度，在許多復雜的水波問題中都得到了很好的驗證和應用[14-17]，而非靜壓模型一般以N-S方程或者RANS方程作為控制方程，因此將VARANS方程與非靜壓模型相結合是非常自然的。綜上，本文基于VARANS方程建立了非靜壓水波模型，可用于模擬島礁地形大粗糙度或類孔隙結構對波浪傳播運動的影響。通過與相應物理模型實驗的對比，以證實在數值模擬中將地形粗糙單元/結構等效為孔隙介質進行計算這一方法的有效性，并研究分析島礁地形表面大粗糙度（孔隙介質）對礁坪上波高和平均水位變化的影響。

2 基于VARANS方程的非靜壓模型

2.1 控制方程

非靜壓模型的一大特點是將自由表面處理成單值函數，從而可以在水面和水底建立貼體網格體系，大大簡化自由表面的捕捉難度并提高計算效率。其中一種方法是將原本建立在直角坐標系（x*,y*,z*）中的控制方程轉化到σ坐標系（x,y,σ）下，相應的σ=1，0分別表示自由表面和水底。將這一轉換過程用于VARANS方程，在σ坐標系下的守恒形式控制方程可表達為

式中，n表示孔隙率，定義為孔隙介質內部一定體積范圍內空腔的占比，當n=1時，式（1）和式（2）即退化為RANS方程；D為總水深，為靜水深（h）和波面高程（η）之和，即D=h+η；U=(Du,Dv,Dw)T，為守恒變量；F、G、H為通量項；Sh表示與水底地形變化相關的底坡源項；Sp表示由水面波動引發的動壓力項；Sν表示湍流耗散項，使用體積平均后的k-ε模型計算湍流[13]；u,v,為σ坐標系下的水質點速度；且與直角坐標系下的垂向速度w存在下述關系：

上述各項的具體表達形式可以參考文獻[10，18]。式（2）中Sr即為描述孔隙介質對流體運動影響的相關項，來源于拓展的Forchheimer關系式，用于封閉體積平均后的RANS動量方程，具體表達式為上式右側每一方向的分量中，第一項為線性項，表示層流狀態下由多孔介質施加的摩擦力；第二項為速度的二次方項，表示湍流運動產生的摩擦力；第三項來源于表示慣性力作用（附加質量）的時間導數項，用以描述流動的非定常狀態。各項對應的系數如下所示：

式中，D50表示多孔介質的中值粒徑；γ是和附加質量相關的經驗系數，通常取0.34；Keulegan-Carpenter系數表示流體水質點運動的特征長度與孔隙介質特征長度的比值[11]，即，其中T為波浪周期。另外，參考Ma等[10]的推薦取值，本文計算取α=200，β=1.1。

與其他采用貼體網格系統的非靜壓模型不同，為了保持邊界條件的連續性，本文并不忽略在自由表面和水底處的切向應力，參考Derakhti等[19]的工作，當水面無風力作用時，對切向方向的速度邊界條件作出如下修正：

2.2 數值格式

采用有限體積和有限差分混合格式離散方程，并參照He等[18]建立的非靜壓水波數值模型，采用三階中心加權無振蕩格式來對單元網格邊界面上的守恒變量進行重構，并使用MUSTA格式計算相應的數值通量，從而建立了Godunov類型的高精度間斷捕捉格式以處理波浪破碎問題。時間步進上采用二階Runge-Kutta方法，并結合顯式的預估-校正格式計算動壓力以修正速度場。

3 數值計算設置與網格收斂性分析

3.1 數值計算地形參數

本文的數值計算參考Buckley等[20]在粗糙（孔隙介質）底床上開展的物理模型實驗進行設置。物理模型實驗的水槽長度為55.0 m，對于數值模擬，為了避免波浪在造波位置出現二次反射，采用域內造波的方式生成波浪，所以相比于模型實驗，水槽在左端增加5.0 m以方便設置造波源和海綿阻尼層（圖1）。造波源距離礁前斜坡坡腳25.5 m，礁前斜坡坡度為1∶5，后接礁坪段長14.0 m，距地面高程為0.7 m，礁坪段末端設置有1∶12的斜坡。圖1中黃黑色虛線段標注的區域表示設置孔隙材料的區域，在相應的物理模型實驗中，使用了邊長為1.8 cm的混泥土立方塊來模擬自然條件下礁坪上的粗糙度，這種方法常用于在物理模型實驗中模擬底部粗糙度[21-22]、植被等[23-24]，也可以用于模擬孔隙材料。在本文的數值計算中，參考物理實驗設置進行換算，在礁前斜坡和礁坪上方相同的區域，設置了孔隙率為0.87，中值粒徑為1.8 cm，高度為2.0 cm的孔隙層材料。

另外，根據模型實驗設置，數值計算域沿程設置了17個波面高程測點（波面以下為負值），各測點坐標位置從左到右依次為-15.0 m、-5.3 m、-3.0 m、-2.0 m、-1.0 m、-0.5 m、0 m、0.34 m、1.0 m、1.6 m、3.0 m、4.2 m、5.8 m、7.5 m、9.0 m和11.0 m。

圖1 礁坪地形及計算域布置示意圖Fig. 1 Schematic of the reef flat and the computational domain layout

3.2 波浪和水位參數

為研究在不同波浪和水位條件下，孔隙介質對礁坪上方波高和增減水的影響是否存在某些共性，因此，在數值計算中一共設置了16個不同的波浪條件和水位組合，其中每一個組合又都包括光滑底床和孔隙底床兩種設置，總計32個工況。域內造波使用Bouws等[25]提出的TMA譜生成波面時間序列，表中的波高為均方波高Hrms。

3.3 網格收斂性分析

本文計算的島礁地形屬于二維剖面，為確定水平x方向的網格尺寸和垂向分層數，本小節對相應方向的網格收斂性進行分析。非靜壓模型能夠通過增加垂向分層數來提高色散性能，有研究認為，有限水深下垂向取3～8層即可滿足精度要求，而水平網格間距應不小于波長的1/40，否則會產生明顯的波幅沿程衰減[10]。本節以表1中的組合1工況（譜峰周期對應的波長為5.5 m）為基礎，設置不同的水平網格大小和垂向分層數，計算沿程各測點波面時間歷程數據的差異性。設置兩組網格參數：（1）保持垂向分層數為8層不變，調整x方向空間步長為0.08 m、0.06 m、0.04 m、0.02 m和0.01 m；（2）保持x方向空間步長為0.01 m不變，調整垂向分層數為2、4、6、8、10，本文所有計算中的垂向分層皆為均勻等間距劃分。另外，為了量化不同網格參數下計算結果的差異性，引入Willmott[26]提出的統計模型：

由圖7可知，隨著硫酸鹽干濕循環作用的增多，試件的相對動彈性模量先提高后降低，并在60次硫酸鹽干濕循環作用處于最高值。原因是侵入混凝土內部的硫酸鹽會與混凝土結晶反應生成石膏等。在較少次數的循環作用下，石膏等表現出的填充作用可以增加混凝土的密實度，故增大了動彈性模量；隨后由于石膏等含量的增多，其表現出的膨脹作用增強，加快了試件內部裂縫的發展，破壞了混凝土的孔隙結構，故動彈性模量值隨之降低[15-17]。

式中，Skill代表兩個數據樣本的相似度；a和b表示不同網格參數的波面數據；ηi,a和 ηi,b分別表示a、b網格參數對應的第i個時刻波面高程；表示網格參數b對應的波面時間序列平均值。另外， 0≤Skill≤1，Skill=1表示兩個波面數據完全一致，等于0則表示差異最大。計算不同網格參數下，沿程各測點波面的Skill值，如圖2所示。

表1 數值實驗波浪和水位參數Table 1 Parameters of wave and water level for the numerical experiments

可以看出，對于x方向不同的網格間距，在礁坪地形前方的常水深處波面差異很小，到達礁前斜坡后隨著波浪的淺化變形，波面差異增大，進入礁坪區域后差異最為明顯，Δx為0.08 m和0.06 m時礁坪處的Skill值約為0.95，但隨著網格的加密，差異減小，Δx為0.02 m和0.01 m的Skill值已非常接近1，可以判斷進一步加密網格不會對計算結果造成影響。而對于不同的垂向分層數，也有著相類似的結果，垂向均勻劃分8層和10層，計算結果間的差別非常小，可以認為垂向分層數超過8時已能獲得很好的收斂效果。

根據上述收斂性分析，在數值計算中，計算域沿x方向離散為6 000個網格，即網格間距取 Δx=0.01 m，垂向均勻劃分10層，計算時間步長由Courant-Friedrichs-Lewy (CFL)條件控制。

4 結果與討論

使用本文建立的基于VARANS方程的非靜壓水波模型，對表1中的16組波浪和水位參數分別在光滑和孔隙地形上的傳播進行了模擬。下面以組合4的計算結果為例進行分析，對應的均方波高為0.12 m，譜峰周期為2.26 s，礁坪水深為0.04 m。

圖2 水平方向網格間距（a）和垂向分層數（b）收斂性分析及地形剖面（c）示意圖Fig. 2 Convergence of grid spacing in horizontal (a) and vertical direction (b) and the sketch of terrain profile (c)

圖3中給出光滑底床和孔隙底床下，離岸處、礁前斜坡頂點和礁坪段的波面時間序列對比，即分別對應x為-15.0 m、0 m和6.0 m。顯然，孔隙底床的存在并沒有對波浪在礁坪地形上的傳播形態造成明顯變化，無論底床光滑與否，呈現出來的都是一個非常典型的簡單斜坡-礁坪地形上波浪的傳播變形過程。波浪從外海（造波源）傳播到礁前斜坡上，由于斜坡對波浪的淺化作用，波形變陡，波高增大，最終在礁坪邊緣附近發生破碎，如圖3b所示，此位置處的波面時間歷程線基本位于0水位線以上，表明在該位置處平均水位上升，即波浪增水。圖3c中測點位于礁坪段，礁坪上水深極淺，由于波能的衰減，波浪形態從曲線上看變現為波高下降明顯，波面歷程線變得相對平坦，“起伏”程度變緩，同時波面過程線完全位于0水位線上方，表明礁坪上方出現了明顯的增水，這可能是礁坪上方波浪高頻成分耗散，低頻成分有所增長在波形上的一種體現。

圖3 組合 4測點波面高程歷程曲線Fig. 3 Time series of simulated wave surface elevation for the Case 4

從圖3中可以看出t=100 s以后的波浪趨于穩定，所以選取組合4工況t為100～500 s（圖3中紅色虛線處截斷）時的波面數據分析均方波高和平均水位（物理實驗中使用的數據時間長度為410 s），如圖4所示。通過和物理模型實驗的測量數據對比，無論是光滑床面還是孔隙床面，可以發現本文模型的計算結果與測量值保持了高度的一致，因此可以認為本研究建立的基于VARANS方程的非靜壓模型可以很好地模擬不規則波在典型礁坪地形上的傳播。對于波高變化，圖3所反映出來的變化與圖2中的波浪形態基本是一致的，即離岸處基本保持穩定，經過斜坡淺化增大，在x=0附近發生破碎；與之對應的是平均水位在島礁地形前方基本為0，在礁坪邊緣前方出現平均水位的下降，但隨后平均水位迅速升高，在礁坪上方體現為波浪增水，同時增水幅值沿礁坪向岸基本不變。

由于島礁前方的深水區域不存在孔隙介質，所以波高和平均水位對于兩種底床都基本一致。但對于孔隙床面，礁前斜坡和礁坪上方的波高都比光滑床面時有所下降，這是由于孔隙介質的存在加大了波浪能量的衰減。對于平均水位而言，孔隙床面情況下的最大減水幅值小于光滑床面，但礁坪上方的波浪增水兩個床面則很接近，事實上，從實驗測量數據上看，個別測點的數據甚至是重合的，說明對于光滑和孔隙兩種底床，礁坪上方的增水幅值并不存在明顯差別。

圖4 組合 4光滑和孔隙底床下均方波高（a）、平均水位（b）和地形剖面（c）沿程變化Fig. 4 Mean square wave height (a) , mean water level (b), and terrain profile (c) for the Case 4 at porous and smooth beds

圖5 孔隙底床和光滑底床在不同區域處均方波高的比較Fig. 5 Comparisons of the mean square wave heights at different regions for the porous and smooth beds

在礁坪上方，各實驗組次對應的均方波高與水深的比值（Hrms/h）如圖6所示。從圖中可以看出在孔隙介質存在的情況下，礁坪上方的波高要小于光滑底床。圖6中同時給出了本文模型數值計算結果和物理模型實驗測量數據的對比情況，顯然，數值計算給出的礁坪上方波高值要略大于實驗測量，但二者的擬合直線（最小二乘法擬合）斜率是接近的。對于數值計算結果，孔隙底床條件下礁坪上方的波高水深比大約為0.4，光滑底床的平均值約為0.3；相應的實驗測量數據則分別為0.36和0.27。

相比于波高的變化，孔隙介質對于平均水位的影響則有所不同。如圖4b所示，各工況在礁坪前方出現最大減水，而光滑底床的減水幅值（絕對值）要大于孔隙底床，表現為圖7中的數據點都位于1∶1斜率直線的上方，對于所有工況，本文模型數值計算結果顯示二者比值的平均值為1.8，物理實驗測量的平均值為1.7。

但是在礁坪上方，平均水位沿空間變化很小，同時平均水位基本一致，無論是數值計算還是實驗測量結果，都反映了這一特點，如圖8所示，數據點基本都分布在1∶1斜率直線上。采用最小二乘法進行擬合，可以發現數值計算結果的擬合直線斜率為0.98，實驗數據對應的斜率為0.94。

圖6 孔隙底床和光滑底床在礁坪上方均方波高和水深比（Hrms/h）Fig. 6 Comparisons of the ratio of mean square wave height to water depth (Hrms/h) on the reef flat for the porous and smooth beds

圖7 孔隙底床和光滑底床條件下最大減水對比Fig. 7 Comparisons of maximum setdown for the porous and smooth beds

Buckley等[20]的物理模型實驗證明，粗糙底床對增水的影響有兩個機制，一是改變了從礁前斜坡到礁坪沿程波浪能量的空間分布，二是同時產生了底部應力。對于前者，孔隙介質帶來的類似底摩擦阻力對波浪速度產生阻滯效果，使得波高和波浪能量出現衰減，相應的導致增水幅值的下降；后者源于孔隙介質（粗糙底床）對波-流場阻力的時間平均，并可用包含速度的二次方項和相應拖曳力系數的阻力定律進行估計[2-3]，這一機制會導致平均水位的上升。因此這兩個機制的影響互相抵消，最終導致礁坪上方的平均水位并未因為孔隙介質的存在而出現顯著變化。

圖8 孔隙底床和光滑底床在礁坪上方波浪增水比較Fig. 8 Comparisons of wave setup on the reef flat for the porous and smooth beds

為更清楚地展現孔隙介質的影響，將孔隙底床各工況均方波高和平均水位相對于光滑底床變化的百分比匯總在表2中。

表2 孔隙底床上均方波高和平均水位相對于光滑底床的變化百分比Table 2 Relative percentage changes of mean square wave heights and average water levels between the porous beds and the smooth beds

一些數值模型通過拖曳力系數（或底摩阻系數）與速度平方的乘積項來計算包含底部粗糙度的流動問題，認為底部粗糙度會造成平均水位上升25%～30%[27-28]，與本文計算結果和相關物理實驗結果之間存在較大差異。對此，有學者指出類似的方法并不能很好地處理此類冠層流動中摩擦阻力對速度的衰減，而相應的一些采用空間平均或沿水深平均的動量方程模型則能獲得很好的效果[24,29]。實際上，VARANS模型也可看作一種冠層流動模型，通過本文的應用，發現其確實對于存在底部摩擦阻力的情況能準確地反映平均水位的變化，相應的，也能反映出孔隙介質的存在確實導致礁坪上方的流速產生了明顯的衰減。如圖9所示，對于光滑底床，礁坪上方靠近水面表現為向岸方向的流動（流速為正），而靠近水底則表現為離岸方向的流動（流速為負），礁坪上方中間的水體則處于兩股反向水流的交界處，流速約為0，Yao等[30]也報導過類似的礁坪上方波生流反向分層現象。孔隙底床的礁坪上方流速出現了明顯的下降，而反向的離岸流動則出現在中間層，另外，靠近水底處，由于位于孔隙介質內部，流速很小，且垂向分布較為均勻。

圖9 組合4光滑底床（a）和孔隙底床（b）礁坪上方水平方向平均流速分布Fig. 9 Distribution of horizontal mean current velocities at the reef flat of Case 4 for the smooth (a) and porous (b) beds

值得指出的是，在目前關于孔隙介質的物理模型實驗中[11,31-32]，孔隙材料的中值粒徑范圍一般為1.5～2.5 cm，孔隙率則集中在0.45～0.52，本文認為，對于實驗中常用的理想均一化的孔隙材料，其孔隙率下限約等于圓球體積與外切立方體的體積比，即大約為0.47。為了進一步研究孔隙介質對礁坪上方增水的影響，下面在組合4原始孔隙率為0.87的基礎上，增加了0.67和0.47兩種孔隙率的計算結果進行對比，如圖10所示。

圖10 組合4不同孔隙率下平均水位對比及地形剖面示意圖Fig. 10 Comparisons of mean water level for Case 4 with varied porosities and the sketch of terrain profile

通過對比可以看出，孔隙率從0.87減小到0.67時，礁坪上方的增水有所上升，但變化不明顯，但當孔隙率進一步減小到0.47時，礁坪上的增水變化更為微小，可以認為達到了一種類似“收斂”的效果。孔隙率變小對應的是孔隙介質的體積比增大，對應在模型實驗中相當于用以模擬孔隙介質或底摩擦阻力的材料的空間布置密度增大，從直觀上來說粗糙程度應當是增大的，但從計算結果看，礁坪上方的增水幅值并無明顯上升，因此可以認為上文提到的底摩擦阻力對波浪增水影響不明顯的論述在此得到進一步的驗證。

5 結論

本文采用基于VARANS方程的非靜壓波浪模型對隨機波浪典型斜坡-礁坪地形上的傳播進行了模擬，重點分析了孔隙介質對礁坪上方均方波高和平均水位變化的影響。通過和相應物理模型實驗結果的對比，證明了本文模型能有效處理波浪和孔隙介質的相互作用，與實驗測量結果保持了很高的一致性。通過分析，發現孔隙介質對波高的衰減作用較為明顯，相比于光滑底床條件，在波浪破碎點附近，孔隙底床的波高下降了12%，礁坪上方波高平均下降28%。對于平均水位，孔隙底床條件下的最大減水幅值減小了43%，同時礁坪上方增水上升6%。相比于波高的變化，礁坪上方增水的變化比較小，但與實驗測量的結果相一致，可以認為孔隙介質（類似于底摩擦阻力的作用）的存在對礁坪上方的增水影響不明顯。另外，孔隙率（0.47～0.87）在一定范圍內變化時，對平均水位的影響不明顯，也從另一個方面證明了孔隙介質對礁坪上方的增水影響很小。