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

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

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

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

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

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

1 數(shù)學模型

1.1 控制方程

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

連續(xù)方程：

(1)

N-S方程：

(2)

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

1.2 湍流模型選擇

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

PT+GT+DkT-εT

(3)

(4)

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

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

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

Vfr=2cη(t)W

(5)

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

(6)

1.4 邊界條件設(shè)置

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

圖1 邊界條件

圖2 模擬試驗段布置

1.5 自由表面追蹤方法

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

(7)

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

(8)

1.6 模型求解與驗證

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

圖3 波面驗證對比

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

圖4 流速驗證對比

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

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

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

(9)

2.1 波面三維運動過程

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

圖5 方形排柱周圍波面三維運動

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

2.2 方形排柱周圍平面速度-渦量場特征

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

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

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

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

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

圖7 H/L對kt的影響

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

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

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

圖8 d/L對kt的影響

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

3 結(jié)語

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

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

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