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

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

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

目前，關(guān)于床面切應(yīng)力的研究多采用物理模型試驗(yàn)的手段。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)力模型試驗(yàn)研究，發(fā)現(xiàn)波浪破碎前后切應(yīng)力具有較大差異。劉清君[7]對(duì)島礁地形上礁緣處和礁坪后方波浪穩(wěn)定區(qū)床面切應(yīng)力進(jìn)行了研究。黃海龍[8]研究了完全直接測量二維波浪和水流共線作用時(shí)床面切應(yīng)力儀的設(shè)計(jì)和制作。

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

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

1 數(shù)學(xué)模型

1.1 控制方程

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

(1)

(2)

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

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

(3)

(4)

(5)

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

1.2 造波和消波邊界

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

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

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

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

(6)

(7)

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

1.4 數(shù)值方法

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

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

2.1 數(shù)值模型建立

假定礁盤地形為梯形平臺(tái)，平臺(tái)高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)對(duì)網(wǎng)格進(jìn)行加密，加密后的網(wǎng)格長度小于波長的1/300，網(wǎng)格高度小于波高的1/20。同時(shí)，為適應(yīng)壁面函數(shù)邊界條件的要求，對(duì)礁坪上近壁面0.01 m范圍進(jìn)行網(wǎng)格加密，加密后網(wǎng)格長度和高度均為0.002 5 m。為保證計(jì)算的穩(wěn)定性和收斂性，計(jì)算時(shí)間步長的選取應(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 計(jì)算工況

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

2.3 模型驗(yàn)證

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

(1)波高驗(yàn)證。

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

(2)斷面流速驗(yàn)證。

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

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中可看出，向岸和離岸的最大平均流速沿水深均存在不同程度的衰減，向岸流的衰減幅度大于離岸流。且隨著波高的增大，向、離岸最大流速分布相對(duì)0值位置的不對(duì)稱性明顯較大，流速分布整體向正方向偏移。數(shù)值模擬研究可以較為準(zhǔn)確地反映出這一現(xiàn)象及沿程不同水深處的流速分布規(guī)律。

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

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

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

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

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

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

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

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

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

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

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

入射波高作為波浪的重要參數(shù)，是波浪在礁坪上變形與破碎的重要影響因素，圖6為△=0.23，周期1.0 s時(shí)，不同入射波高條件下，向岸和離岸切應(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ā)生破碎，破碎點(diǎn)處波高增加，使得向岸切應(yīng)力增加，而離岸切應(yīng)力主要由波谷作用時(shí)離岸回流引起，波谷傳播至水深突變的礁緣處水體回落，會(huì)在礁緣處產(chǎn)生劇烈的離岸流。

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

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

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

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

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

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

綜上可見，周期對(duì)穩(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ù)值波浪水槽，對(duì)礁盤地形上波浪傳播進(jìn)行數(shù)值模擬，分析了相對(duì)水深、入射波高和波周期對(duì)礁坪上的波浪底部切應(yīng)力沿程變化的影響，得到主要結(jié)論如下：

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

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

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

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