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

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

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

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

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

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

1 數(shù)學模型

1.1 控制方程

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

(1)

(2)

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

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

(3)

(4)

(5)

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

1.2 造波和消波邊界

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

1.3 固壁邊界及底部切應(yīng)力

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

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

(6)

(7)

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

1.4 數(shù)值方法

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

2 數(shù)值模型建立與驗證

2.1 數(shù)值模型建立

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

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

圖1 數(shù)值域網(wǎng)格及邊界條件設(shè)置Fig.1 Numerical grides and boundary condition of the numerical domain

2.2 計算工況

數(shù)值計算中考慮波周期、波高和礁坪上相對水深對波浪底部切應(yīng)力的影響，試驗波浪采用規(guī)則波，波高選取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，不同參數(shù)之間進行組合，共計10組計算工況。

2.3 模型驗證

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

(1)波高驗證。

參考模型試驗中波高測點位置，在數(shù)值模型中分別選取距離坡腳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值位置的不對稱性明顯較大，流速分布整體向正方向偏移。數(shù)值模擬研究可以較為準確地反映出這一現(xiàn)象及沿程不同水深處的流速分布規(guī)律。

(3)底部切應(yīng)力驗證。

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

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 數(shù)值計算結(jié)果分析

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

3.1 礁坪上相水深對底部切應(yīng)力影響

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

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

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

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

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

3.2 入射波高對底部切應(yīng)力影響

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

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

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

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

3.3 入射波周期對底部切應(yīng)力影響

圖7為不同波周期條件下，向岸和離岸切應(yīng)力沿程變化關(guān)系圖。

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

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

綜上可見，周期對穩(wěn)定后的向岸切應(yīng)力影響大于離岸切應(yīng)力。

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

4 結(jié)論

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

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

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

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

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