基于OpenFOAM的礁坪上波浪底部切應力數值模擬研究

黃 哲，李巖汀，胡乾明，劉清君，孫天霆，王登婷*

(1.南京水利科學研究院，南京 210024；2.水文水資源與水利工程科學國家重點實驗室，南京 210029；3.河海大學 港口海岸與近海工程學院，南京 210098)

礁盤地形通常由陡峭的礁前斜坡和水深較淺的礁坪組成，波浪在礁緣或礁坪上傳播會發生劇烈的變形和破碎。波浪破碎區是泥沙起動及物質輸運的重要區域，因此，開展該區域床面切應力研究，分析底部切應力隨水深、入射波高和波周期等因素變化規律，對于研究礁盤地形上的生態環境和地形演變具有重要的科學和現實意義。

目前，關于床面切應力的研究多采用物理模型試驗的手段。Mirfenderesk[1]使用應力板分別測量了光滑和粗糙床面的波浪床面切應力。Seelen[2]使用應力板測量了孤立波作用下的床面切應力。Arnskov[3]使用熱膜探針測量了波、流共同作用下的床面切應力。Tanaka[4]應用嵌平式熱膜探針在波浪水槽內測量了橢圓余弦作用下層流邊界層內床面切應力。徐華等[5]采用柔性熱膜式壁面切應力傳感器測量了波、流共同作用下床面切應力。郝思禹[6]在坡度1：5的斜坡上開展了切應力模型試驗研究，發現波浪破碎前后切應力具有較大差異。劉清君[7]對島礁地形上礁緣處和礁坪后方波浪穩定區床面切應力進行了研究。黃海龍[8]研究了完全直接測量二維波浪和水流共線作用時床面切應力儀的設計和制作。

關于床面切應力的數值模擬研究一般通過計算區域流場或紊動場，然后根據流速與切應力之間關系計算得到底部切應力，常見計算方法有二次阻力律法、對數擬合流速剖面法等。Jonsson[9]提出波浪摩阻系數定義式，通過邊界層外近底最大自由流速計算床面切應力。Cox[10]在此基礎上引入近底床面瞬時流速，建立了瞬時床面切應力二次阻力律公式。張慶河[11]通過對大量波浪邊界層和床面切應力數據進行分析，建立了可用于波浪作用下層流和紊流邊界層床面摩阻系數計算式。流速剖面法假定近底流速沿水深近似滿足對數分布，通過近底區域內流速擬合得到摩阻流速，在通過摩阻流速計算床面切應力[12]。Yuan等[13]在振蕩水槽中對比各種計算床面切應力計算方法，認為流速剖面法準確性最好。

結合國內外研究情況，以往波浪底部切應力的研究大都集中于平底地形或單一的斜坡地形，針對復合地形的試驗測點布置較少，不足以分析波浪底部切應力在礁盤地形上的沿程變化規律[7]。為此，本文基于OpenFOAM開源代碼中兩相流求解器OlaFlow建立礁盤地形的數值波浪水槽，對礁盤地形上的波浪傳播進行數值模擬，分析水深、入射波高和波周期對礁坪上波浪底部切應力沿程變化的影響。

1 數學模型

1.1 控制方程

OpenFOAM中OlaFlow兩相流求解器采用的控制方程為三維VARANS方程

(1)

(2)

式中：ui為xi方向的速度分量；ρ為流體密度；t為時間；p*為動壓力；g為重力加速度；Xj為控制體中心點的位置分量；μtot為有效動力粘度系數，μtot=μ+μt，μ為粘度系數，μt為紊動粘度系數，取決于所選用的紊流模型。

采用SSTk-ω紊流模型[14]確定VARANS方程中的紊動粘度系數μt，具體如下

(3)

(4)

(5)

式中：k為紊動動能；ω為紊動能耗散率；S為應變率張量；F1、F2為混合函數；紊流模型中各參數取值如下：β*=0.09，α1=5/9，β1=3/40，α2=0.44，β2=0.082 8。

1.2 造波和消波邊界

在計算域入口使用速度入口邊界條件[15]，速度入口邊界條件通過設定波面和流體速度進進行造波。根據LeMehaute波浪理論適用范圍[16]，結合本文所模擬的波浪要素，采用二階斯托克斯波理論進行造波。在計算域入口和出口處設置主動消波邊界[17]進行消波，消波邊界通過產生與入射波速度相同方向相反的速度進行消波。主動消波方法相比于松弛域消波或多孔介質消波等被動消波方法，可以顯著減小計算域大小，節省計算資源。

1.3 固壁邊界及底部切應力

數值水槽中一般將固壁邊界設為不可滑移邊界條件，但這樣規定邊界條件就必須將微分方程對黏性底層積分。由于黏性底層中的流速梯度很大，為了得到理想的數值解，必須在近壁面區域布置極細密的網格，極大的增加了計算成本，也容易造成網格畸形導致計算發散。為了解決上述問題，模型中引入壁面函數，建立了近壁面網格點至壁面之間的流速非線性變化模型。底部切應力可以利用切應力的定義式并結合壁面函數求出。

壁面切應力按式(6)計算，式中的系數vnew由式(7)確定。

(6)

(7)

式中：uc為靠近壁面第一層網格中心點平行于壁面的流速；uw為壁面上的流體速度，可取0；κ為卡門常數，取0.41；E為常數，取9.8；y+為靠近壁面第一層網格中心點到壁面的無量綱垂直距離。

1.4 數值方法

OpenFOAM使用有限體積法進行離散求解。本文各控制方程中，時間導數項采用隱式歐拉格式，梯度項采用高斯線性插值，對流項采用Gauss limited Linear V1格式，拉普拉斯項采用高斯線性修正格式，速度壓力耦合方程采用PIMPLE(PISO-SIMPLE)算法求解。PIMPLE算法是將單個時間步長內流動看作是穩態流動使用SIMPLES算法進行求解，在時間步進上采用PISO算法進行求解。

2 數值模型建立與驗證

2.1 數值模型建立

假定礁盤地形為梯形平臺，平臺高0.5 m，迎浪斜坡坡度為1:1，右側礁坪長度為8 m，礁盤總長8.5 m。

礁盤地形數值波浪水槽網格及邊界條件示意圖見圖1。水槽長14 m，寬1 m，高1 m，礁盤坡腳設置在距離造波邊界3倍波長處。為能夠精確捕捉自由水面和避免波浪沿程衰減，在靜水面上下3倍波高范圍內對網格進行加密，加密后的網格長度小于波長的1/300，網格高度小于波高的1/20。同時，為適應壁面函數邊界條件的要求，對礁坪上近壁面0.01 m范圍進行網格加密，加密后網格長度和高度均為0.002 5 m。為保證計算的穩定性和收斂性，計算時間步長的選取應使得Courant數小于1，本文取波周期的1/1 000。

圖1 數值域網格及邊界條件設置Fig.1 Numerical grides and boundary condition of the numerical domain

2.2 計算工況

數值計算中考慮波周期、波高和礁坪上相對水深對波浪底部切應力的影響，試驗波浪采用規則波，波高選取0.05 m、0.06 m、0.07 m和0.08 m，波周期選取1.0 s、1.5 s、2.0 s和2.5 s，礁前水深選取0.6 m、0.65 m、0.7 m和0.75 m，不同參數之間進行組合，共計10組計算工況。

2.3 模型驗證

利用礁盤地形上波浪傳播物模試驗結果[7]，通過波面過程線、最大平均流速沿水深分布和底部切應力等方面對比，驗證所建礁盤地形數值波浪水槽模型。

(1)波高驗證。

參考模型試驗中波高測點位置，在數值模型中分別選取距離坡腳0.5 m、1.0 m、2.0 m、6.5 m的4個測點，給出波浪歷時曲線對比圖，如圖2所示。圖中可以看出，波高過程計算值與實測值吻合良好，表明該模型能夠有效地模擬波浪在礁盤地形上的淺水變形和波浪破碎。

(2)斷面流速驗證。

為了進一步驗證模型的合理性，分別對斷面最大平均流速沿水深分布進行了驗證。其中斷面最大平均流速選取距離坡腳1.3 m和2.3 m兩個斷面，測點水深分別為-0.2h2、-0.4h2、-0.6h2、-0.73h2和-0.95h2，流速驗證情況如圖3所示，其中向岸為正，離岸為負。

X=0.5 m X=1.0 m

X=2.0 m X=6.5 m2-a H=0.05 m

X=0.5 m X=1.0 m

X=2.0 m X=6.5 m2-b H=0.07 m

X=1.3 m X=2.3 mX=1.3 mX=2.3 m3-a H=0.05 m3-b H=0.07 m

圖3中可看出，向岸和離岸的最大平均流速沿水深均存在不同程度的衰減，向岸流的衰減幅度大于離岸流。且隨著波高的增大，向、離岸最大流速分布相對0值位置的不對稱性明顯較大，流速分布整體向正方向偏移。數值模擬研究可以較為準確地反映出這一現象及沿程不同水深處的流速分布規律。

(3)底部切應力驗證。

底部切應力驗證選取的兩個測點分別位于距離坡腳0.6 m和4.0 m處，圖4給出了剪切力隨時間的變化對比圖，圖中向岸切應力為正，離岸切應力為負，切應力計算值與實測值吻合較好，極值較為接近。

X=0.6 m X=4.0 m4-a H=0.05 m

X=0.6 m X=4.0 m4-b H=0.07 m

3 數值計算結果分析

為了分析波浪底部切應力在礁坪上的沿程變化，從礁緣處開始，每隔0.5 m設置數據采集點，用于采集流速和紊動粘性系數，然后通過計算得到底部切應力。

3.1 礁坪上相水深對底部切應力影響

礁坪上相對水深△(△=h2/h1)是影響礁坪上波浪破碎和底部切應力的重要因素，其中h2為礁坪上水深，h1為礁前水深。

圖5反映了入射波要素為H=0.06 m，T=1.0 s時，不同△條件下向岸和離岸切應力沿程變化情況。整體來看，向岸和離岸切應力均隨著△的減小而呈現出增大趨勢，主要原因是隨著△的減小，波浪傳播受礁盤地形影響增大，波浪發生淺水變形，波峰變得尖陡而波谷變得平坦，非線性效應逐漸增強，當波高達到破碎值后波浪發生破碎，產生的紊動水體會擴散并影響近底流速，進而引起切應力增大。

圖5中向岸和離岸切應力在礁坪上的沿程變化情況可以看出，當△較大時，向岸和離岸切應力在礁坪上沿程變化幅度相對較小，而隨著△的減小，切應力沿程變化幅度增大。結合△=0.17時切應力沿程變化情況可見，向岸和離岸切應力均隨與礁緣間距離的增大呈現先增大后減小的趨勢，且在礁緣處的4倍波長左右(X=1.5 m)達到峰值，主要是由于波浪傳播至該處時發生卷破，水舌投入水中后引起水體發生劇烈紊動，導致離岸和向岸切應力增大。

波浪在礁坪上部傳播一段距離后趨于穩定，除△=0.23，其他△工況下，穩定后的切應力大小相近，表明礁坪上水深對波浪穩定后的底部切應力影響相對較小，其值一般由入射波浪條件控制。

5-a 向岸切應力 5-b 離岸切應力圖5 沿程切應力隨相對水深變化(H=0.06 m，T=1.0 s)

3.2 入射波高對底部切應力影響

入射波高作為波浪的重要參數，是波浪在礁坪上變形與破碎的重要影響因素，圖6為△=0.23，周期1.0 s時，不同入射波高條件下，向岸和離岸切應力沿程變化圖。

總體上看，向岸和離岸切應力均隨著入射波高的增加而增大。向岸切應力最大值大都出現在X=1.0 m附近，而離岸切應力主要出現在礁緣處，且離岸切應力的最大值約為向岸切應力的2.5倍至3倍。這是由于向岸切應力主要受波高影響，波浪傳播至礁坪上發生破碎，破碎點處波高增加，使得向岸切應力增加，而離岸切應力主要由波谷作用時離岸回流引起，波谷傳播至水深突變的礁緣處水體回落，會在礁緣處產生劇烈的離岸流。

6-a 向岸切應力 6-b 離岸切應力圖6 沿程切應力隨入射波高變化(T=1.0 s，h2=0.15 m)

圖6-a中，不同入射波高對波浪穩定后切應力有一定影響，入射波高越大，向岸切應力越大，而圖6-b中，不同波高穩定后離岸切應力差別相對較小，入射波高對穩定后離岸切應力幾乎無影響。

3.3 入射波周期對底部切應力影響

圖7為不同波周期條件下，向岸和離岸切應力沿程變化關系圖。

圖7-a可見向岸切應力最大值出現在礁緣后方破波點附近。隨著入射波周期的增加，礁坪上最大切應力增大且最大值所在位置向礁緣后方移動。在X>7.5 m之后，向岸切應力沿程增加，這是因為波浪破碎后會形成段波向后方傳播，行進波的形成影響近底流速，使向岸切應力增加。圖中可以看出波周期越大，穩定后的向岸切應力也越大。

由圖7-b可見，離岸切應力在礁緣處達到最大值后迅速衰減，不同于向岸切應力，離岸切應力在礁坪上保持沿程衰減，在波浪傳播達到穩定后，大周期波浪相應的離岸切應力略大，但不同周期下波浪穩定處的切應力差值整體較小。

綜上可見，周期對穩定后的向岸切應力影響大于離岸切應力。

7-a 向岸切應力 7-b 離岸切應力圖7 沿程切應力隨入射波周期變化(H=0.06 m，h2=0.15 m)

4 結論

本文基于OpenFOAM中兩相流求解器OlaFlow建立礁盤地形數值波浪水槽，對礁盤地形上波浪傳播進行數值模擬，分析了相對水深、入射波高和波周期對礁坪上的波浪底部切應力沿程變化的影響，得到主要結論如下：

(1)礁坪上向岸和離岸切應力隨相對水深的減小而增大，隨入射波高和入射波周期的增大而增大；

(2)礁坪上波浪穩定后的向岸底部切應力隨入射波高和入射波周期的增大而增大，受相對水深影響較小。

(3)礁坪上波浪穩定后的離岸底部切應力受入射波高、入射波周期和相對水深影響均較小。

(4)向岸切應力最大值主要出現在礁緣后方，隨入射波周期增大，向岸切應力最大值出現位置向礁緣后方移動；離岸切應力最大值主要出現在礁緣處。