劉 文，賀昌海，惠建偉，劉 全

(1.武漢大學(xué)水資源與水電工程科學(xué)國家重點實驗室，武漢 430072；2.中鐵二十局集團有限公司，西安 710016)

數(shù)值模擬是繼物理模型試驗之后研究施工導(dǎo)流中水力學(xué)問題的又一重要手段，目前已有較多關(guān)于導(dǎo)流流場及局部范圍內(nèi)的數(shù)值模擬研究。楊研等[1]基于非結(jié)構(gòu)ELCIRC模型對航道水流條件進行數(shù)值模擬，分析了通航條件下的局部礙航流態(tài)；戚藍等[2]對基于精確河道地形的溢流壩泄流進行三維數(shù)值模擬,數(shù)值模擬結(jié)果與試驗值吻合良好,驗證了數(shù)值模擬的可靠性；羅永欽等[3]結(jié)合大比尺單體模型試驗成果，應(yīng)用分段計算方法對泄洪洞摻氣減蝕問題進行了三維紊流數(shù)值模擬分析。驗證了泄洪洞各道摻氣坎布置合理且有穩(wěn)定摻氣空腔，摻氣坎體型基本合理；王月華等[4]基于Flow-3D軟件同時結(jié)合水工模型試驗具體分析了消能池的水力特性，比較全面地反映了消力池的水流和消能情況；程香菊[5]等用水汽兩相流混合模型模擬計算了階梯溢流壩面水氣兩相流的流速、摻氣濃度、初始摻氣點和壓力等特征參數(shù)，實測結(jié)果驗證了數(shù)值模擬的可靠性、合理性；高學(xué)平等[6]基于雙流體歐拉法對溢洪道摻氣挑坎摻氣水流進行水氣兩相流三維數(shù)值模擬研究，研究了溢洪道沿程及斷面的摻氣濃度分布規(guī)律等特性，驗證了雙流體歐拉法研究溢洪道摻氣水流的可行性。

前人研究證明數(shù)值模擬和模型試驗結(jié)果吻合較好。但是，這些成果要么只是針對簡化的整體流場，要么只是考慮單個水工建筑物的數(shù)值模擬，沒有包括考慮水流摻氣的導(dǎo)流泄水建筑物的整體導(dǎo)流流場的模擬成果。

本文以某分期導(dǎo)流工程為例，建立三維導(dǎo)流模型，考慮溢洪道(作為導(dǎo)流泄水建筑物)消力池內(nèi)的水流摻氣影響，與模型試驗結(jié)果對比分析,以期獲得更加合理、精確的數(shù)值模擬成果。

1 三維模型

某工程由左右岸土堤、左岸心墻壩、河床心墻壩、溢洪道、灌溉取水口及廠房組成。左右岸土堤Hmax=24.8 m，Lleft=940 m，Lright=4 735.84 m；左岸心墻壩Hmax=44.8 m，L=190.5 m；河床心墻壩Hmax=51.3 m，L=608.5 m；溢洪道最大高度54.8 m，寬141 m，設(shè)6個泄洪底孔、2個泄洪表孔，消力池段長112.99 m，消力坎段長28 m。溢洪道底孔溢流面高程485.00 m，孔口尺寸8.0 m(寬)×8.5 m(高)。表孔溢流面高程507.00 m，孔寬15 m。工程采用兩期導(dǎo)流法施工，二期導(dǎo)流標準為100 a一遇，流量8 800 m3/s，上游組合圍堰設(shè)計頂高程505.00 m，在二期導(dǎo)流階段，溢洪道表孔溢流堰體不施工，預(yù)留缺口作為導(dǎo)流泄水建筑物，缺口底部高程485.00 m(見圖1)。

圖1 導(dǎo)流平面布置Fig.1 Layout of diversion project

1.1 地形模型

按1∶1的比尺，采用CATIA 軟件生成地形三維模型[7,9](見圖2)。

圖2 工程三維河道地形Fig.2 The 3D river terrain

1.2 溢洪道及建筑物模型

進入CATIA 創(chuàng)成式外形設(shè)計工作平臺，建立溢洪道的開挖面，在零件設(shè)計工作平臺通過建立一系列開挖面進行開挖，形成溢洪道[8]，在開挖后的三維模型上建立圍堰和泄水建筑物(比尺1∶1)(見圖3)。

圖3 導(dǎo)流工程三維整體模型Fig.3 The 3D model of diversion project

消力池及閘室細部結(jié)構(gòu)見圖4、圖5。

圖4 消力池細部Fig.4 The 3D model of stilling pool

圖5 閘室進口Fig.5 The inlet of the sluice chamber

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

2.1 紊流模型

本文研究的導(dǎo)流問題地形起伏多變，消力池區(qū)域水流流態(tài)復(fù)雜、流線彎曲程度大，故采用RNGk-ε模型[9]?？刂品匠倘缦?。

連續(xù)方程：

(1)

動量方程：

(2)

紊動能k方程：

(3)

紊動能耗散率ε方程：

(4)

式中：ui為流速分量；Ai、gi、fi分別為三維直角坐標方向上可流動的面積分數(shù)、重力加速度和黏滯力；VF是可流動的體積分數(shù)；ρ是流體密度；p是作用在流體微元上的壓力；μt是紊動黏滯系數(shù)；Gk是紊動能k的產(chǎn)生項：σk、σε是紊動能和耗散率對應(yīng)的Prandtl數(shù)，均為1.39；經(jīng)驗常數(shù)Cε1、Cε2分別取1.42、1.68。

2.2 卷氣模型

消力池中伴隨著強烈的卷氣，采用Flow-3D中的卷氣模型(Air Entrainment Model)?？紤]到流體體積膨脹及卷入氣體的浮力效應(yīng)，卷氣模型中還要使用密度變化方程。卷氣模型基于3種因素：由紊流產(chǎn)生的擾動、重力及表面張力，紊流是卷氣過程產(chǎn)生的主要因素。流體表面摻氣是因為紊流渦體在流動過程中將表面空氣卷入流體中，其主要取決于紊流強度是否克服重力和表面張力組成的表面穩(wěn)定力[10]。

紊流渦體的特征尺寸:

(5)

單位體積的紊動能：

Pt=ρQ

(6)

表面穩(wěn)定力:

(7)

單位時間卷氣體積量:

(8)

式中：Q為紊動能；D為耗散函數(shù)；在RNG紊流模型中，cnu的值為0.085，用來描述表面擾動的特征；ρ是液體密度；gn是重力加速度在自由表面法線方向的分量；Lt為流體單元升高高度；σ是表面張力系數(shù)；As是表面面積；Cair是一個經(jīng)驗參數(shù),表示單位時間內(nèi)有一部分表面積卷入空氣，初步假設(shè)值為0.5，文獻中已驗證了其在數(shù)值模擬中的合理性[10,11]；當紊動能Pt小于表面穩(wěn)定力Pd時，δV值為0。

2.3 漂移通量模型

在漂移通量模型(Drift-flux Model)中,兩相速度以某一混合速度為基礎(chǔ),通過定義漂移速度分別得到氣、液各相速度，在下面的守恒方程中將相對速度視為漂移速度，通過附加的氣相連續(xù)性方程來描述氣液兩相流動[12]。在本文中，漂移通量模型模擬在摻氣過程中氣泡相對于水流的漂移運動以及空氣因為浮力效應(yīng)而離開流體的現(xiàn)象。

混合相連續(xù)性方程：

(9)

氣相連續(xù)性方程：

(10)

動量方程：

(11)

(12)

3 三維數(shù)值模擬

為了保證計算效率、穩(wěn)定性和準確性，采用均勻的立方體網(wǎng)格，其縱橫比為1。為更精確地模擬消力池內(nèi)的流動，消力池處采用較小的網(wǎng)格單元，模型均采用結(jié)構(gòu)化網(wǎng)格。網(wǎng)格劃分包含3個網(wǎng)格塊(見圖6)，覆蓋范圍為下游0-800 m至上游0+1 000 m處，消力池范圍內(nèi)(網(wǎng)格塊②)立方體網(wǎng)格單元邊長1.2 m，上下游立方體網(wǎng)格單元邊長1.8 m,網(wǎng)格總數(shù)量約2 000 萬個。

圖6 網(wǎng)格劃分Fig.6 Mesh generation

模型的上游邊界設(shè)置為流量邊界，給定上游水位和流量，下游邊界設(shè)置為壓力邊界，給定相應(yīng)的下游水位。所有網(wǎng)格塊底部設(shè)為壁面無滑移邊界，其他邊界保留默認設(shè)置，為對稱邊界，水面設(shè)置為自由液面。

初始條件設(shè)置為各個工況的初始上游水位，壓力為靜水壓。壓力求解選擇隱式GMRES算法，同時采用RNG紊流模型、卷氣模型、漂移通量模型、變密度模型(Density Evaluation Model)，并用FLOW-3D基于結(jié)構(gòu)化矩形網(wǎng)格的 FAVOR 方法模擬實體及其獨有的真實3步Tru-VOF 方法追蹤自由表面。計算工況見表1。

表1 數(shù)值模擬工況Tab.1 Condition of numeric simulation

4 數(shù)值模擬與試驗成果對比分析

4.1 流 態(tài)

流態(tài)對比結(jié)果見圖7、圖8，可看出數(shù)值模擬和模型試驗?zāi)M的缺口進口水流都比較平靜。工況2，試驗的閘室出口流速較大，流態(tài)紊亂，出口處有明顯的水躍產(chǎn)生，水流摻氣明顯，數(shù)值模擬流態(tài)與之相符合。

圖7 消力池內(nèi)水流流態(tài)對比(工況1)Fig.7 Contrast of flow regime (condition 1)

圖8 消力池內(nèi)水流流態(tài)對比(工況2)Fig.8 Contrast of flow regime (condition 2)

4.2 水 位

溢洪道測點水位對比參見表2。從表2中可以看出，試驗水位與計算水位總體吻合較好，但工況1與工況2樁號0-142 m處水位相差較大。其原因可能是，消力池在樁號0-50.8 m至0-165.6 m之間, 測點0-142 m處于消力池尾部，一方面消力池內(nèi)流態(tài)紊亂，水流紊動導(dǎo)致消力池尾部水位波動較大，測量時同一點水位不同時刻差異明顯，且空氣被卷入水流，基本不存在界限分明的自由水面；另一方面，數(shù)值模擬計算達到穩(wěn)定狀態(tài)后，不同時刻該點的水位也能得到不同的值。

表2 溢洪道測點水位對比 m

計算值比試驗值稍大，以工況2為例，計算值的誤差平均值0.49，誤差最大值0.98。上下游圍堰前水位波動變幅很小，上游臨時圍堰(設(shè)計頂高程505.00 m)和下游圍堰(設(shè)計頂高程496.50 m)具有較大的安全裕度，能夠安全擋水。

在數(shù)值模擬中分別計算考慮摻氣影響和未考慮摻氣影響的工況。在工況1中，考慮摻氣影響時(加入卷氣模塊)的上游圍堰迎水側(cè)水位為492.28 m，未考慮摻氣影響的上游圍堰迎水側(cè)水位為492.24 m。在工況2中，考慮摻氣影響的上游圍堰迎水側(cè)水位為504.27 m，未考慮摻氣影響的上游圍堰迎水側(cè)水位也為504.27 m，均低于上游圍堰設(shè)計高程505.00 m。由此可知，在本工程中考慮溢洪道消力池內(nèi)的水流摻氣影響與否，對施工導(dǎo)流的上游圍堰高度基本沒有影響。

4.3 流 速

溢洪道中心線及圍堰附近測點流速對比見圖9，樁號0+300、0+203、0-182、0-300是指圍堰附近測點，其余樁號則表示溢洪道中心線各測點。圍堰附近測點計算流速值與模型試驗值結(jié)果幾乎一致。工況1和工況2中斷面0-99.6 m和斷面0-160 m處出現(xiàn)了流速誤差值大于1 m/s的情況，工況1中計算值的相對誤差平均值28.69%，最大誤差-1.43 m/s。工況2中計算值的相對誤差平均值25.43%，最大誤差1.31 m/s。究其原因，這2個斷面正好位于消力池內(nèi)水流強烈摻氣范圍內(nèi)，消力池內(nèi)流態(tài)復(fù)雜，流速脈動性強。

圖9 溢洪道中心線及圍堰附近測點流速結(jié)果對比Fig.9 Contrast of flow velocity at measurement points of the spillway center line & cofferdam

4.4 壓強分布

溢洪道中心線測點壓強對比見圖10。上下游壓強分布較均勻，無劇烈變化，閘室出口和消力池內(nèi)壓力變化較劇烈且消力池底部壓力較大。

圖10 溢洪道中心線測點壓強對比Fig.10 Contrast of pressure at measurement points along the spillway center line

工況1中計算值的相對誤差平均值13%，工況2中計算值的相對誤差平均值7%?？傮w而言，考慮水流摻氣時,大部分壓強結(jié)果吻合較好，但消力池中的壓強有一定偏差。

4.5 水流摻氣率

整體流場中水流摻氣濃度分布見圖11，據(jù)此可確定水流摻氣范圍。

圖11 水流摻氣濃度平面分布Fig.11 Plane figure of air entrainment concentration distribution

由水流摻氣濃度截面分布圖12可知，斷面摻氣濃度自消力池底部向水面方向逐漸增加，與實際相符。且隨著流量增大，初始摻氣點后移，確定初始摻氣點的位置，有助于確定非摻氣區(qū)的范圍，并采取相應(yīng)的措施避免可能發(fā)生的空蝕破壞[13]。

圖12 水流摻氣濃度截面分布Fig.12 Profile figure of air entrainment concentration distribution

由數(shù)值模擬結(jié)果，求得7個測點處的表面摻氣率(見表3)。

表3 表面摻氣率Tab.3 Air entrainment ratio on the surface

5 結(jié) 論

本文結(jié)合FAVOR實體邊界模擬技術(shù)模擬復(fù)雜地形，基于FLOW-3D軟件的Tru-VOF技術(shù)追蹤水流自由表面，綜合紊流模型、卷氣模型和漂移通量模型,對某實際分期導(dǎo)流工程全流場進行了三維數(shù)值模擬，得到了較好的結(jié)果，驗證了上游圍堰設(shè)計高程的合理性，加入卷氣模塊模擬三維流場得到了較好的流速分布、壓強分布、水流摻氣濃度和摻氣范圍，為類似導(dǎo)流工程的流場計算提供了一種理論和方法。