基于OpenFOAM的礁坪上波浪底部切應(yīng)力數(shù)值模擬研究

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

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

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

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

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

結(jié)合國(guó)內(nèi)外研究情況，以往波浪底部切應(yīng)力的研究大都集中于平底地形或單一的斜坡地形，針對(duì)復(fù)合地形的試驗(yàn)測(cè)點(diǎn)布置較少，不足以分析波浪底部切應(yīng)力在礁盤(pán)地形上的沿程變化規(guī)律[7]。為此，本文基于OpenFOAM開(kāi)源代碼中兩相流求解器OlaFlow建立礁盤(pán)地形的數(shù)值波浪水槽，對(duì)礁盤(pán)地形上的波浪傳播進(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]，速度入口邊界條件通過(guò)設(shè)定波面和流體速度進(jìn)進(jìn)行造波。根據(jù)LeMehaute波浪理論適用范圍[16]，結(jié)合本文所模擬的波浪要素，采用二階斯托克斯波理論進(jìn)行造波。在計(jì)算域入口和出口處設(shè)置主動(dòng)消波邊界[17]進(jìn)行消波，消波邊界通過(guò)產(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ā)散。為了解決上述問(wèn)題，模型中引入壁面函數(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；κ為卡門(mén)常數(shù)，取0.41；E為常數(shù)，取9.8；y+為靠近壁面第一層網(wǎng)格中心點(diǎn)到壁面的無(wú)量綱垂直距離。

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í)間步長(zhǎng)內(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ù)值模型建立

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

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

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

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

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

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

為了進(jìn)一步驗(yàn)證模型的合理性，分別對(duì)斷面最大平均流速沿水深分布進(jìn)行了驗(yàn)證。其中斷面最大平均流速選取距離坡腳1.3 m和2.3 m兩個(gè)斷面，測(cè)點(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ì)稱(chēng)性明顯較大，流速分布整體向正方向偏移。數(shù)值模擬研究可以較為準(zhǔn)確地反映出這一現(xiàn)象及沿程不同水深處的流速分布規(guī)律。

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

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

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)力在礁坪上的沿程變化，從礁緣處開(kāi)始，每隔0.5 m設(shè)置數(shù)據(jù)采集點(diǎn)，用于采集流速和紊動(dòng)粘性系數(shù)，然后通過(guò)計(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)力沿程變化情況。整體來(lái)看，向岸和離岸切應(yīng)力均隨著△的減小而呈現(xiàn)出增大趨勢(shì)，主要原因是隨著△的減小，波浪傳播受礁盤(pá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)力沿程變化情況可見(jiàn)，向岸和離岸切應(yīng)力均隨與礁緣間距離的增大呈現(xiàn)先增大后減小的趨勢(shì)，且在礁緣處的4倍波長(zhǎng)左右(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)力幾乎無(wú)影響。

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

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

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

綜上可見(jiàn)，周期對(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建立礁盤(pán)地形數(shù)值波浪水槽，對(duì)礁盤(pán)地形上波浪傳播進(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)在礁緣處。