規(guī)則波作用下方形排柱附近水流特性數(shù)值模擬*

王澤能，劉世建，呂文舒，陳 瑞

(1.廣西北部灣投資集團(tuán)有限公司，廣西 南寧530022；2.廣西交科集團(tuán)有限公司，廣西 南寧 530022；3.北京工業(yè)大學(xué)，北京 100124；4.長(zhǎng)沙理工大學(xué) 土木工程學(xué)院，湖南 長(zhǎng)沙 410015)

波浪與方柱相互作用是流體力學(xué)領(lǐng)域的經(jīng)典問(wèn)題，相關(guān)研究開(kāi)展較早。萬(wàn)德成等[1]運(yùn)用VOF方法追蹤自由表面，較早地建立了孤立波翻越直立方柱流動(dòng)的數(shù)值水槽；李玉成等[2]通過(guò)試驗(yàn)研究，對(duì)方柱在規(guī)則波作用下的流場(chǎng)特征進(jìn)行了初步探討。

近年來(lái)，隨著多種透空式防波堤[3-4]、透水丁壩[5]、透水潛壩等[6]新型建筑物的提出，波浪、水流作用下透空結(jié)構(gòu)周圍水動(dòng)力特性正引起相關(guān)研究學(xué)者與工程技術(shù)人員的注意。結(jié)合恰當(dāng)?shù)耐牧髂Ｐ烷_(kāi)展數(shù)值模擬研究，可充分了解此類結(jié)構(gòu)周圍流場(chǎng)特征組成部分。曹洪建等[7]運(yùn)用naoe-FOAM-SJTU求解器，模擬三維潰壩波在方柱周圍劇烈流動(dòng)的過(guò)程，計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)吻合較好。蔣昌波等[8]、姚宇等[9]先后在孤立波條件下對(duì)圓形排柱周圍水動(dòng)力特性進(jìn)行了數(shù)值模擬研究。海顯盛等[10]采用Boltzmann方法成功模擬了彎曲河道方形橋墩周圍的紊動(dòng)特征。陳浩民等[11]運(yùn)用有限元方法離散Berkhoff緩坡方程計(jì)算非淹沒(méi)雙直立圓柱周圍波高分布。殷銘簡(jiǎn)等[12]采用OpenFoam建立數(shù)值波浪水槽，探討孤立波與密排樁防波堤相互作用的數(shù)值模擬研究。從近期相關(guān)成果可以看出，排柱式防波堤[13]、海洋鉆井平臺(tái)、跨海大橋與海上風(fēng)電場(chǎng)等[14]都有排柱類結(jié)構(gòu)的運(yùn)用。然而前人大多是研究單個(gè)方柱在孤立波背景下的流動(dòng)特性，或是圓形排柱周圍流動(dòng)特性，對(duì)于規(guī)則波作用下方形排柱周圍的水流特性認(rèn)識(shí)尚有不足。關(guān)于規(guī)則波同方形排柱的相互作用問(wèn)題尚有待進(jìn)一步研究。

本文基于FLOW-3D計(jì)算流體力學(xué)軟件，采用RNGk-ε方程控制的N-S方程配合VOF自由表面追蹤方法建立波浪-方形排柱相互作用水槽，對(duì)規(guī)則波作用下方形排柱附近的流動(dòng)特征進(jìn)行數(shù)值模擬研究，可為相關(guān)實(shí)際工程提供參考。

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

1.1 控制方程

采用連續(xù)方程和雷諾時(shí)均N-S方程作為波浪作用下方形排柱附近周圍流體運(yùn)動(dòng)的控制方程。

連續(xù)方程：

(1)

N-S方程：

(2)

式中：t為時(shí)間；ρ為水的密度；P為壓力；u、v、w分別為各方向上的速度分量；Ax、Ay、Az分別為各方向上的面積分?jǐn)?shù)；Gx、Gy、Gz分別為各方向上的重力加速度；VF為體積分?jǐn)?shù)；fx、fy、fz分別為各方向上的黏滯力加速度。

1.2 湍流模型選擇

大渦模擬和RNGk-ε湍流模型均具有準(zhǔn)確模擬結(jié)構(gòu)物周圍流場(chǎng)特征的能力，但大渦模擬需要占用大量的計(jì)算資源，并需要非常精細(xì)的網(wǎng)格，綜合計(jì)算能力與處理成本高。選用RNGk-ε湍流模型，模型紊動(dòng)動(dòng)能kT和紊動(dòng)耗散率εT的表達(dá)式如下：

PT+GT+DkT-εT

(3)

(4)

式中：PT為速度梯度產(chǎn)生的紊動(dòng)動(dòng)能；GT為浮力產(chǎn)生的紊動(dòng)動(dòng)能；C1、C3分別取1.42和0.20，C2由kT和PT計(jì)算得到；DkT和Dε分別為對(duì)應(yīng)面積分?jǐn)?shù)和體積分?jǐn)?shù)上的紊動(dòng)擴(kuò)散項(xiàng)。

1.3 數(shù)值水槽源項(xiàng)設(shè)定

參考Lin[15]提出的質(zhì)量源造波法，在數(shù)值水槽中添加一個(gè)長(zhǎng)L、寬W、高H的質(zhì)量源模型，通過(guò)給定質(zhì)量源模型上隨時(shí)間變化的體積流率來(lái)實(shí)現(xiàn)波浪模擬，具體的體積流率Vfr表達(dá)式如下：

Vfr=2cη(t)W

(5)

式中：c為波浪傳播速度；W為水槽寬度；η(t)為波面方程。在波浪模擬的前3個(gè)周期內(nèi)，對(duì)體積流率疊加包絡(luò)函數(shù)R，預(yù)防計(jì)算發(fā)散。

(6)

1.4 邊界條件設(shè)置

在模擬波浪傳播時(shí)，底部為無(wú)滑移固壁邊界；頂部為對(duì)稱邊界；兩側(cè)為對(duì)稱邊界；為了保持水槽內(nèi)的水體質(zhì)量守恒，水槽出入口均設(shè)置成壓力邊界，并給定靜水壓強(qiáng)和靜水位高程值。質(zhì)量源模型兩側(cè)均有波形生成，在水槽出入口設(shè)置1～2倍波長(zhǎng)海綿層區(qū)域。試驗(yàn)段布置在水槽中段x=10 m處，并將加密嵌套網(wǎng)格設(shè)為對(duì)稱邊界，數(shù)值水槽布置見(jiàn)圖1，方形排柱模擬段幾何參數(shù)見(jiàn)圖2。

圖1 邊界條件

圖2 模擬試驗(yàn)段布置

1.5 自由表面追蹤方法

FLOW-3D軟件采用三步流體體積法(Tru-VOF)來(lái)追蹤自由液面的流動(dòng)情況。當(dāng)界面參數(shù)aw=0，當(dāng)前網(wǎng)格中視為不存在液體；當(dāng)界面參數(shù)aw的取值在0～1，表示當(dāng)前網(wǎng)格氣液共存；當(dāng)界面參數(shù)aw=1，表示當(dāng)前網(wǎng)格被液體充滿，不存在氣體。氣、液交界面追蹤由以下方程來(lái)求解：

(7)

VOF法流體體積函數(shù)為F=F(x，y，z，t)，其表達(dá)形式如下：

(8)

1.6 模型求解與驗(yàn)證

采用有限差分法對(duì)控制方程進(jìn)行離散處理，以GMERS方法配合自適應(yīng)步長(zhǎng)方法進(jìn)行求解。將質(zhì)量源波面方程設(shè)置為水深d=0.2 m、波高H=0.01 m、周期T=1 s的規(guī)則波方程進(jìn)行數(shù)值模擬。數(shù)值水槽全長(zhǎng)20 m，將質(zhì)量源模型設(shè)置在距數(shù)值水槽左邊界x=5 m處。在水槽中段設(shè)置大小為0.01 m的均勻網(wǎng)格，在數(shù)值水槽兩端設(shè)置3 m的海綿層并配合設(shè)置變步長(zhǎng)漸變網(wǎng)格段波吸收波浪，共計(jì)約550萬(wàn)個(gè)網(wǎng)格，數(shù)值試驗(yàn)均在這一驗(yàn)證網(wǎng)格下開(kāi)展。提取水槽中段波面數(shù)據(jù)同數(shù)值解進(jìn)行對(duì)比。圖3顯示，同理論值相比，所建立的波浪數(shù)值水槽波面情況具有良好的重復(fù)性。

圖3 波面驗(yàn)證對(duì)比

為驗(yàn)證數(shù)值水槽模擬結(jié)構(gòu)周圍流速過(guò)程的準(zhǔn)確性，參考王珍等[16]進(jìn)行的規(guī)則波作用下墩柱周圍流動(dòng)特性試驗(yàn)，在數(shù)值水槽中建立相同圓形墩柱結(jié)構(gòu)，設(shè)置水深d=0.6 m、波高H=0.17 m、周期T=2.0 s，選取水深d=6.8 cm處墩柱正面特征點(diǎn)流速進(jìn)行比較。圖4顯示，數(shù)值水槽計(jì)算得到的水平流速與橫向流速同實(shí)測(cè)值相比誤差在合理范圍內(nèi)，證明所建立的波浪數(shù)值水槽具有良好的流速計(jì)算能力。

圖4 流速驗(yàn)證對(duì)比

2 方形排柱周圍波浪傳播特性數(shù)值模擬結(jié)果分析

設(shè)定多種波高、周期、水深開(kāi)展多種工況數(shù)值模擬。通過(guò)波浪自由面特征與流速、渦量等定量指標(biāo)對(duì)規(guī)則波作用下方形排柱周圍波浪傳播特性研究。

方形排柱設(shè)置足夠高度，在模擬過(guò)程中不發(fā)生頂面越浪。在方形排柱后x=11 m處設(shè)置浪高監(jiān)測(cè)點(diǎn)記錄排柱周圍水面線歷時(shí)特征，透射系數(shù)是評(píng)價(jià)結(jié)構(gòu)后方波浪傳播特性的主要指標(biāo)，根據(jù)透射波高Ht和入射波高Hi計(jì)算透射系數(shù)Kt，見(jiàn)式(9)。

(9)

2.1 波面三維運(yùn)動(dòng)過(guò)程

選取典型波況(d=0.5 m，T=1.2 s，H=0.14 m)，根據(jù)數(shù)值水槽嵌套加密區(qū)的模擬結(jié)果，在波面穩(wěn)定后，繪制方形排柱周圍三維波面運(yùn)動(dòng)圖，在特征時(shí)刻對(duì)方形排柱周圍一個(gè)波周期T內(nèi)的波面運(yùn)動(dòng)情況進(jìn)行觀測(cè)，見(jiàn)圖5。波浪沿x軸正向方形排柱布置區(qū)傳播。在圖5a)、b)所示前半周期內(nèi)，波浪逐漸進(jìn)入方形排柱影響區(qū)域范圍內(nèi)，由于方形排柱的阻滯作用，方形排柱前端的反射波面初步形成。

圖5 方形排柱周圍波面三維運(yùn)動(dòng)

圖5c)～e)所示波浪周期內(nèi)，波峰到達(dá)排柱前端，入射波波峰同排柱前反射波疊加在排柱前形成壅高，當(dāng)入射波波峰運(yùn)動(dòng)至排柱軸線時(shí)，入射波波面在排柱的作用下發(fā)生分離，水體從排柱間的間隙區(qū)域涌入排柱后方。排柱后方水體相互摻混形成復(fù)雜的排柱后方波面分布特征。當(dāng)波峰線越過(guò)排柱軸線后，波面逐漸恢復(fù)，由于結(jié)構(gòu)阻滯與水體摻混的能量耗散，波高降低，以透射波形式繼續(xù)向后傳播，繼而波谷區(qū)域通過(guò)排柱軸線區(qū)域，在波浪負(fù)向速度梯度的作用下，水體回流，壅水出現(xiàn)在樁后，見(jiàn)圖5f)。從波面模擬情況上看，所建立的數(shù)值水槽三維表面追蹤能力強(qiáng)，可以精細(xì)刻畫(huà)排柱結(jié)構(gòu)周圍流場(chǎng)特征。

2.2 方形排柱周圍平面速度-渦量場(chǎng)特征

方形排柱周圍波面的多樣變化，提示排柱周圍水體發(fā)生了復(fù)雜的波浪-結(jié)構(gòu)相互作用。同波浪與單體結(jié)構(gòu)物作用不同，方形排柱的幾何特征決定其周圍流場(chǎng)特性有其自身特征。選擇典型波況(d=0.5 m，T=1.2 s，H=0.14 m)，復(fù)合提取相對(duì)水深z/d=0.6處xy平面上一個(gè)波周期T內(nèi)典型時(shí)刻的渦量與速度矢量分布，對(duì)方柱周圍流場(chǎng)特性進(jìn)行分析，見(jiàn)圖6。

圖6 方形排柱周圍水平流速渦量

數(shù)值模擬采用線性波方程控制，在波浪靜水面附近，水質(zhì)點(diǎn)水平運(yùn)動(dòng)速度為0，在波峰和波谷處取得水質(zhì)點(diǎn)水平運(yùn)動(dòng)最大速度。如圖6a)～c)所示，當(dāng)波谷經(jīng)過(guò)排柱軸線分布區(qū)域時(shí)，排柱間負(fù)向流速不斷增大，且由于排柱壓縮波面，使得排柱間距中的水體在波谷經(jīng)過(guò)排柱軸線時(shí)獲得最大負(fù)向流速；同時(shí)，渦量峰值出現(xiàn)在方柱前端角隅處。隨著波浪推進(jìn)，負(fù)向流速梯度減小，當(dāng)波峰半周期經(jīng)過(guò)時(shí)渦量峰值在方柱角隅處反向，隨波浪水體的傳播出現(xiàn)在方柱后端角隅處。在一個(gè)波浪周期內(nèi)，最大渦量的分布位置在方柱角隅處交替出現(xiàn)，充分反映了波浪的周期性運(yùn)動(dòng)特征與方形排柱的結(jié)構(gòu)特性。

2.3 波陡H/L對(duì)方形排柱波浪透射系數(shù)Kt的影響

取固定水深d=0.5 m，在不同波周期T=1.0～1.4 s和不同波高H=0.08～0.14 m條件下開(kāi)展數(shù)值模擬研究，分析波陡H/L對(duì)方形排柱后方透射系數(shù)Kt的影響，見(jiàn)圖7。

圖7 H/L對(duì)kt的影響

在數(shù)值模擬試驗(yàn)范圍內(nèi)，3種相對(duì)水深(d/L=0.194、0.244、0.330)條件下的透射系數(shù)Kt隨波陡H/L的變化略有波動(dòng)，但整體評(píng)價(jià)影響并不顯著。在固定水深條件下，隨著相對(duì)水深的增大，波長(zhǎng)減小則意味著波周期減小，較大波長(zhǎng)的波浪經(jīng)過(guò)方形排柱時(shí)其透射系數(shù)較大。

2.4 相對(duì)水深d/L對(duì)方形排柱波浪透射系數(shù)Kt的影響

取固定波高H=0.10 m，在不同波周期T=1.0～1.6 s和不同水深d=0.45～0.55 m條件下開(kāi)展數(shù)值模擬研究，分析相對(duì)水深d/L對(duì)方形排柱后方透射系數(shù)Kt的影響，見(jiàn)圖8。

圖8 d/L對(duì)kt的影響

在數(shù)值模擬試驗(yàn)范圍內(nèi)，3種相對(duì)波高(H/d=0.182、0.200、0.222)條件下的透射系數(shù)Kt隨著相對(duì)水深d/L的增大逐步減小。在固定波高條件下，隨著水深增加，透射系數(shù)Kt亦有所增大。對(duì)于類似方形排柱的設(shè)計(jì)計(jì)算，應(yīng)更關(guān)注此類工程的工作水深情況。

3 結(jié)語(yǔ)

1)與相關(guān)試驗(yàn)數(shù)據(jù)對(duì)比，基于RNGk-ε湍流模型封閉與VOF方法追蹤自由表面所建立的三維質(zhì)量源造波數(shù)值水槽工作性能良好。

2)利用建立的數(shù)值水槽對(duì)方形排柱周圍波面演化情況進(jìn)行數(shù)值模擬研究，排柱周圍波面變化表現(xiàn)出較強(qiáng)的三維特征。水平速度分布與渦量場(chǎng)的演化特征分析，揭示了規(guī)則波作用下方形排柱周圍流場(chǎng)呈周期性變化規(guī)律，最大渦量值在方形排柱單柱的上下游角隅區(qū)交替出現(xiàn)。

3)利用無(wú)量綱參數(shù)波陡H/L與相對(duì)水深d/L分析方形排柱后方波浪透射系數(shù)Kt的變化規(guī)律，在數(shù)值模擬試驗(yàn)范圍內(nèi)，透射系數(shù)Kt隨波陡H/L的變化略有波動(dòng)，但整體評(píng)價(jià)影響并不顯著；透射系數(shù)Kt隨相對(duì)水深d/L的增大逐步減小。