考慮水流摻氣的施工導(dǎo)流工程流場(chǎng)數(shù)值模擬

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

(1.武漢大學(xué)水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430072；2.中鐵二十局集團(tuán)有限公司，西安 710016)

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

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

本文以某分期導(dǎo)流工程為例，建立三維導(dǎo)流模型，考慮溢洪道(作為導(dǎo)流泄水建筑物)消力池內(nèi)的水流摻氣影響，與模型試驗(yàn)結(jié)果對(duì)比分析,以期獲得更加合理、精確的數(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個(gè)泄洪底孔、2個(gè)泄洪表孔，消力池段長(zhǎng)112.99 m，消力坎段長(zhǎng)28 m。溢洪道底孔溢流面高程485.00 m，孔口尺寸8.0 m(寬)×8.5 m(高)。表孔溢流面高程507.00 m，孔寬15 m。工程采用兩期導(dǎo)流法施工，二期導(dǎo)流標(biāo)準(zhǔn)為100 a一遇，流量8 800 m3/s，上游組合圍堰設(shè)計(jì)頂高程505.00 m，在二期導(dǎo)流階段，溢洪道表孔溢流堰體不施工，預(yù)留缺口作為導(dǎo)流泄水建筑物，缺口底部高程485.00 m(見(jiàn)圖1)。

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

1.1 地形模型

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

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

1.2 溢洪道及建筑物模型

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

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

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

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

圖5 閘室進(jìn)口Fig.5 The inlet of the sluice chamber

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

2.1 紊流模型

本文研究的導(dǎo)流問(wèn)題地形起伏多變，消力池區(qū)域水流流態(tài)復(fù)雜、流線彎曲程度大，故采用RNGk-ε模型[9]。控制方程如下。

連續(xù)方程：

(1)

動(dòng)量方程：

(2)

紊動(dòng)能k方程：

(3)

紊動(dòng)能耗散率ε方程：

(4)

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

2.2 卷氣模型

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

紊流渦體的特征尺寸:

(5)

單位體積的紊動(dòng)能：

Pt=ρQ

(6)

表面穩(wěn)定力:

(7)

單位時(shí)間卷氣體積量:

(8)

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

2.3 漂移通量模型

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

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

(9)

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

(10)

動(dòng)量方程：

(11)

(12)

3 三維數(shù)值模擬

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

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

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

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

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

4 數(shù)值模擬與試驗(yàn)成果對(duì)比分析

4.1 流 態(tài)

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

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

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

4.2 水 位

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

表2 溢洪道測(cè)點(diǎn)水位對(duì)比 m

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

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

4.3 流 速

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

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

4.4 壓強(qiáng)分布

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

圖10 溢洪道中心線測(cè)點(diǎn)壓強(qiáng)對(duì)比Fig.10 Contrast of pressure at measurement points along the spillway center line

工況1中計(jì)算值的相對(duì)誤差平均值13%，工況2中計(jì)算值的相對(duì)誤差平均值7%。總體而言，考慮水流摻氣時(shí),大部分壓強(qiáng)結(jié)果吻合較好，但消力池中的壓強(qiáng)有一定偏差。

4.5 水流摻氣率

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

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

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

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

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

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

5 結(jié) 論

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