練繼建， 任盼紅， 劉東明， 何軍齡

(1.天津大學(xué) 水利工程仿真與安全國家重點實驗室， 天津 300350； 2.天津大學(xué) 建筑工程學(xué)院， 天津 300350)

1 研究背景

水體中摻入空氣形成水氣二相流是水利工程中常見的現(xiàn)象。水體摻氣對泄水建筑物有利有弊，有利的方面主要有：水體摻氣可以降低脈動壓強[1]，還可以減免空蝕破壞，通過采用摻氣設(shè)施使過水建筑物表面水流發(fā)生強迫摻氣從而保護泄水建筑物[2-4]。不利的方面主要有：水體摻氣會使水深增加[5-6]，若是對水體摻氣量估計不足，會導(dǎo)致開敞式溢洪道邊墻設(shè)計高度偏低，從而使得摻氣水流越過邊墻，危害泄水建筑物安全，或?qū)е乱绾槎春托购槎炊错斢喾鶞p小，通氣不暢，從而造成泄水建筑物的損壞[7-8]。泄洪洞與溢洪道的不同之處在于泄洪洞中水體摻入的空氣并非直接來自于大氣，而是大氣先經(jīng)補氣洞流入泄洪洞中，隨后部分空氣摻入水體中，余下部分則隨水流流出泄洪洞[6,8]。影響洞內(nèi)摻氣特性的因素不僅包括來流流速、來流水深及水流流態(tài)等水力因素，還包括泄洪洞幾何形態(tài)和閘門型式等結(jié)構(gòu)布置因素[9]。

多位學(xué)者基于模型試驗對洞內(nèi)摻氣進行了研究。Hohermuth等[10-12]基于來流水頭達(dá)30.0 m、水槽長度達(dá)20.6 m的模型試驗對洞內(nèi)摻氣進行了研究，發(fā)現(xiàn)空氣需求量隨著閘門后收縮斷面弗勞德數(shù)的增加而增加，且摻氣濃度與水流流速等水力參數(shù)均表現(xiàn)出與明渠流摻氣相似的特性。Aydin[13]通過建立閘門關(guān)閉后水流的非恒定數(shù)學(xué)模型，得到了通風(fēng)井內(nèi)的空氣流量以及閘門后壓力隨時間的變化規(guī)律。Speerli等[14]基于隧道物理模型試驗，研究了不同工況下隧道內(nèi)高速水流的水-氣特征，并分析了水體中空氣濃度的發(fā)展過程，得到了計算最大空氣濃度的公式。岳書波等[15]對高速明流泄洪洞進行了分析與研究，發(fā)現(xiàn)閘室處的進氣量隨著閘門后水流弗勞德數(shù)的增大而增大。何佳等[16]對龍落尾泄洪洞起始段摻氣設(shè)施進行了試驗研究，發(fā)現(xiàn)隨著工作水頭的增加，過坎弗勞德數(shù)先減小后增大。物理模型試驗一般僅滿足重力相似準(zhǔn)則，難以反映原型水體的表面張力和粘滯力，造成模型試驗存在難以忽略的縮尺效應(yīng)[17-18]。

近年來，計算機技術(shù)快速發(fā)展，數(shù)值模擬使得研究摻氣水流的手段更加豐富[19]。黃宗柳[20]和Tang[21]采用數(shù)值模擬方法對明渠自摻氣水流與破碎波摻氣進行了計算，得到的摻氣濃度與試驗數(shù)據(jù)[22]吻合良好。針對隧道內(nèi)不同條件下氣流量的計算，Yazdi等[23]基于FLUENT對平板閘門后隧洞中的高速水流進行了三維數(shù)值模擬研究，并推導(dǎo)了關(guān)于隧洞內(nèi)補氣量和水流上方空氣速度的公式。張法星等[24]采用標(biāo)準(zhǔn)k-ε紊流模型和VOF(volume of fluent model)方法對糯扎渡水電站泄洪洞進行了模擬計算，驗證了k-ε紊流模型結(jié)合VOF方法可以用于預(yù)測大型泄洪洞的需氣量。鄧軍等[25]采用FLUENT對某泄洪洞進行了計算，得到的壓強和通氣設(shè)施進氣量等參數(shù)的計算結(jié)果和模型試驗對比較好。此外，李玲等[26]采用k-ε紊流模型和VOF方法對溪洛渡泄洪洞進行了數(shù)值模擬研究，得到了沿程水面線以及底板中心線的壓力分布。羅永欽等[27]對溪洛渡水電站的一條泄洪洞進行了數(shù)值模擬研究，結(jié)合模型試驗數(shù)據(jù)，認(rèn)為分段計算方法對于高流速、大梯度泄洪洞問題來說可以提高計算效率。Wei等[28]對隧道中的自由面水流供氣特性進行了數(shù)值模擬研究，在不考慮摻氣的條件下，提出了水流上方平均風(fēng)速的計算公式，其平均誤差為±25%。王孝群[8]采用改進的拖曳力模型對補氣洞的補氣量進行了計算，總體誤差近20%。劉嘉夫等[29]與李明達(dá)等[30]分別基于FLOW-3D、FLUENT對摻氣水流進行了數(shù)值模擬研究，計算結(jié)果同樣與試驗數(shù)據(jù)存在誤差。針對以往數(shù)值模擬結(jié)果存在誤差這一問題，本文基于自主開發(fā)的NEWTANK[31]數(shù)值模型，構(gòu)建了研究摻氣問題的混合-濃度模型，通過有限差分方法進行離散并求解，對洞內(nèi)摻氣試驗進行了數(shù)值模擬計算，并探究了不同來流水頭對水流流速和摻氣濃度等參數(shù)的影響。

2 數(shù)值模型建立

2.1 混合模型

本研究針對的是洞內(nèi)摻氣水流，水體和氣體的速度均小于聲速，因此采用不可壓縮流體的控制方程。采用的混合模型首先對除重力項外的其他項進行了水-氣兩相之間的平均，獲得混合流場和混合壓力場；然后采用基于雷諾平均的方法簡化流場內(nèi)的紊動發(fā)展過程，這種方法計算量小，運算速度快?；旌夏Ｐ途唧w表達(dá)式如下：

連續(xù)性方程：

(1)

動量方程：

(2)

式中：下標(biāo)i、j=1，2，3分別表示x，y，z3個方向；um為混合流體的速度，m/s；t為時間，s；pm為混合流體的壓強，Pa；ρ0為流體背景密度，kg/m3；ν為分子的運動黏滯系數(shù)，m2/s；νt為混合流體的紊動黏滯系數(shù)，通過紊流模型來獲取，m2/s；ρm為混合流體的密度，kg/m3；g為重力加速度，m/s2。

重力項可通過下式進行化簡，體現(xiàn)了摻氣水流中空氣與水體間的相互作用：

(3)

式中：C為每個網(wǎng)格中空氣所占體積比，即摻氣濃度。當(dāng)C=1時，表示該網(wǎng)格中充滿了空氣；當(dāng)C=0時，表示該網(wǎng)格中充滿了水體；當(dāng)0

將公式(3)代入公式(2)中，動量方程最終可表示為如下形式：

(4)

為方便表示，后續(xù)水-氣混合相關(guān)參數(shù)均無下標(biāo)“m”，但依舊代表各相流體的混合量。

紊流模型采用應(yīng)用較為廣泛的k-ε兩方程模型，具體表達(dá)如下：

(5)

(6)

(7)

(8)

式中：k為紊動動能，m2/s2；ε為紊動耗散率，m2/s3；Pb為浮力制造項，m2/s3；Prt取常數(shù)為0.85；σk、σε、C1ε、C2ε、C3ε、Cd均為經(jīng)驗系數(shù)，分別取值為σk=1.0、σε=1.0、C1ε=1.44、C2ε=1.92、C3ε=1.44、Cd=0.09。

2.2 濃度模型

本文采用的濃度模型[32]不考慮氣體在水體中的溶解，亦不考慮各相流體間的化學(xué)反應(yīng)，用于追蹤自由液面位置以及獲取摻氣水流內(nèi)部的摻氣濃度分布信息，具體表達(dá)式如下：

(9)

式中：σa為氣體的擴散系數(shù)，表示紊動和氣泡間的相互作用，在此其值取為1；ub為各相流體間的相對滑移速度，m/s。

Haberman等[33]提出了單個氣泡在靜水中的上浮速度計算公式：

(10)

式中：ub0為氣泡在靜水中的上浮速度，m/s；Da為氣泡直徑，m；σt為表面張力，在此式中其值取為0.073 2 N/m。

然而，氣泡在高速水流中的上浮速度的影響因素要比靜水中復(fù)雜很多。當(dāng)水流紊動作用劇烈時，紊動會將大氣泡撕裂成小氣泡，氣泡互相接觸、干擾，導(dǎo)致氣泡間的碰撞、變形、破碎與聚并，造成氣泡動能的損失，從而造成氣泡上浮速度減小。在浮力作用下氣泡傾向于向水體自由面上浮，在湍流擴散作用下氣泡傾向于從高濃度區(qū)域移動到低濃度區(qū)域，當(dāng)周圍水體中均含有空氣時，摻氣濃度梯度減小，造成氣泡上浮速度減小。為盡量體現(xiàn)紊動水流中的氣泡上浮情況，本文建立的數(shù)值模型采用的氣泡上浮速度同時考慮了紊動作用和摻氣濃度的影響[20]，計算公式如下所示：

(11)

2.3 數(shù)值模型計算流程

本數(shù)值模型的計算流程如圖1所示。

圖1 數(shù)值模型計算流程圖

3 模型設(shè)置與驗證

3.1 模型設(shè)置

高水頭結(jié)構(gòu)物出口處常配置閘門，從閘門后流出的水流會帶走大量空氣，易導(dǎo)致負(fù)壓，從而引發(fā)空蝕空化等問題。通過通風(fēng)口提供充足的空氣可有效緩解該問題，然而，盡管通風(fēng)口對閘下出流的安全很重要，但關(guān)于這方面的設(shè)計研究卻很少。Hohermuth[12]在蘇黎世聯(lián)邦理工學(xué)院水力學(xué)、水文學(xué)和冰川實驗室進行了關(guān)于明流隧洞內(nèi)水流摻氣特性的物理模型試驗研究，試驗中，進口水深通過不帶閘門槽的矩形尖頂閘門進行控制，閘門最大高度為0.25 m，明流隧洞寬度為0.2 m，高度為0.3 m，隧洞坡度為0.04。靠近閘門處連接一個圓形通風(fēng)口，其直徑為0.1 m。本文對上游來流水頭為10 m(相應(yīng)的單寬流量為0.88 m2/s)，隧洞長度為6.6 m，閘門開度為40%的試驗工況進行了數(shù)值模擬計算。

依據(jù)試驗中物理模型的設(shè)置以及水流的來流條件，將數(shù)值算例計算域設(shè)置為6.8 m×0.5 m，并采用均勻網(wǎng)格進行離散。在x和z方向分別設(shè)置均勻網(wǎng)格170和100個，最小網(wǎng)格尺寸分別為0.040和0.005 m。在數(shù)值模擬過程中，采用坐標(biāo)變換的方法，將原來傾斜的隧洞變成水平。z方向位于閘門處，垂直于槽底，x方向與隧洞底部重合。上游通風(fēng)口設(shè)置在距離進口0.2 m的位置，由于通風(fēng)口采用圓形截面，而隧洞采用矩形截面，因此可通過圓形面積除以隧洞寬度從而求得通風(fēng)口的等效長度，即0.04 m。進口依據(jù)工況條件給定水深和水流流速，其中流速大小通過閘門開度和流量計算求得；出口處給定自由出流；底部采用墻邊界條件；頂部給定壓力進口。洞內(nèi)水流摻氣數(shù)值模型邊界條件示意圖見圖2。

圖2 洞內(nèi)水流摻氣數(shù)值模型邊界條件示意圖

3.2 摻氣濃度驗證

數(shù)值模擬得到的摻氣濃度和流場分布如圖3所示。由圖3可見，水流自閘門流出后，由于水體流速較高，紊動作用劇烈，空氣逐漸摻入水體中，使水深逐漸增大。隨著水流流向下游，摻氣量逐漸增多，水深繼續(xù)增大，最后達(dá)到相對穩(wěn)定狀態(tài)。

圖3 摻氣濃度和流場分布數(shù)值模擬結(jié)果

物理模型中，在x=1.00 m、x=2.00 m、x=2.99 m、x=3.99 m、x=4.98 m和x=5.98 m 6個斷面進行摻氣濃度和水流流速的測量。摻氣濃度的數(shù)值模擬結(jié)果與試驗數(shù)據(jù)對比如圖4所示。

圖4 摻氣濃度數(shù)值模擬結(jié)果與試驗數(shù)據(jù)對比

由圖4可看出，摻氣濃度的數(shù)值模擬結(jié)果與試驗數(shù)據(jù)總體上較為吻合，僅局部存在差異。靠近進口位置的隧洞底部摻氣濃度與試驗數(shù)據(jù)吻合程度較高，而靠近出口位置處，底部的摻氣濃度數(shù)值模擬結(jié)果大于試驗結(jié)果，分析其原因是由于底部氣泡上浮速度被低估，導(dǎo)致氣泡不能及時逸出水面，造成數(shù)值模擬計算的摻氣濃度比試驗值偏大。

x=1.00 m斷面數(shù)值模擬計算得到的摻氣水深略大于試驗值，其他斷面處數(shù)值模擬計算得到的摻氣水深均略低于試驗值。差異的原因可能來自兩個方面：一是距閘門較遠(yuǎn)位置靠近自由面的氣泡上浮速度在數(shù)值模擬計算中被高估，使得氣泡過早地逸出水面，造成水深偏低；二是相較于數(shù)值模擬計算得到的相對光滑自由表面，在模型試驗實際過程中，由于水流流速較大，紊動程度較高，摻氣過程的發(fā)生使得自由面不規(guī)則，波動的自由面摻雜著飛濺而出的小水滴，從而造成摻氣水深偏大。

各斷面處摻氣濃度分布形狀基本一致：在低摻氣濃度區(qū)域(C<0.2),摻氣濃度隨深度的變化較為明顯；在高摻氣濃度區(qū)域(C≥0.2)，即靠近水-氣自由面位置，摻氣濃度變化梯度較小，隨水深變化不明顯。x=2.00 m處靠近隧洞底部的試驗數(shù)據(jù)局部突然增大是由于模型試驗自身突出的接縫而產(chǎn)生的，在數(shù)值模型的設(shè)置中此處無法正確建模，因而造成數(shù)值模擬計算結(jié)果和試驗數(shù)據(jù)存在較大差異。

3.3 水流流速驗證

6個斷面處水流流速的數(shù)值模擬結(jié)果與試驗數(shù)據(jù)對比如圖5所示。

由圖5可看出，在距離入口較近的位置，數(shù)值模擬得到的速度比試驗值偏小；隨著與入口位置距離的增加，差異逐漸減小；直至在靠近出口處(如x=5.98 m)差異基本消失，此時數(shù)值模擬得到的流速和試驗數(shù)據(jù)幾乎重疊。氣泡大小隨時間和空間的變化而變化，而本數(shù)值模型選取的固定氣泡直徑或許更加接近下游水流中的真實情況，因而造成這種差異沿程變化的現(xiàn)象。模型試驗只給出了摻氣水流的混合流速，而未給出空氣的速度，由圖5數(shù)值模擬計算結(jié)果可明顯看出：僅在摻氣水流上方一定高度，空氣是流向下游出口方向的，而在靠近隧洞上壁面位置，存在與水流流向相反的氣流，且靠近上壁面位置的氣流速度隨著與入口位置距離的增大而逐漸增大。在x=1.00 m處空氣速度達(dá)到了-2 m/s左右，而在x=5.98 m處空氣速度接近-5 m/s，隨著與入口位置距離的逐漸增加，摻氣水流的水深逐漸增加，因而逆向補氣區(qū)域高度逐漸降低，由于逆向補氣量是恒定的，因而造成氣流速度逐漸增大。

圖5 水流流速數(shù)值模擬結(jié)果與試驗數(shù)據(jù)對比

前文中圖3顯示了此工況下計算域內(nèi)的流場分布情況?？諝庾陨嫌瓮L(fēng)孔進入隧洞后，靠近水體的部分空氣從通氣孔位置流向下游出口，與此同時還存在隧洞出口處的逆向補氣，這相反方向的兩股氣流在洞內(nèi)形成一個大的逆時針氣流旋渦。進口處流出的水流因流速高、紊動大，洞內(nèi)的部分空氣被摻入水體中，并帶向下游，導(dǎo)致負(fù)壓產(chǎn)生，在壓差的作用下，上游通風(fēng)口將外部大氣帶入洞內(nèi)。若不考慮水流的拖曳作用，在壓差的作用下空氣會從出口位置進入隧洞；若不考慮壓差作用，水流的拖曳作用將使空氣向出口方向運動。該工況下洞內(nèi)余幅較大，摻氣水流的拖曳能力較為有限，僅能抵抗一部分壓差作用，因而在拖曳作用下僅帶走靠近摻氣水體的部分空氣，靠近洞頂位置的空氣將在壓差的作用下由隧洞出口進入洞內(nèi)，形成一個大的氣流循環(huán)，對應(yīng)的氣流循環(huán)示意圖如圖6所示。

圖6 隧洞內(nèi)氣流循環(huán)示意圖

4 來流水頭對水流流速和摻氣濃度的影響

基本的水氣特征參數(shù)包括水流流速和摻氣濃度，設(shè)定同一閘門開度，通過改變來流水頭(對應(yīng)單寬流量見表1)對摻氣濃度、平均摻氣濃度和水流流速等參數(shù)沿水流方向的變化規(guī)律進行了探究。

表1 閘門開度為40%時不同來流水頭對應(yīng)單寬流量

4.1 對水流流速的影響

對于給定的閘門開度，當(dāng)改變來流水頭時各斷面水流流速分布如圖7所示。由圖7可以看出，當(dāng)閘門開度不變，來流水頭逐漸增大時，水流流速的分布形式相似；當(dāng)閘門開度一定，來流水深保持一致時隨著來流水頭的增加，來流流量逐漸增大，水流流速也隨之增大。

圖7 閘門開度為40%時不同來流水頭He下各斷面水流流速分布

4.2 對摻氣濃度的影響

保持閘門開度為40%時，不同來流水頭各斷面摻氣濃度分布如圖8所示。由圖8可以看出，對于給定來流水頭，距離閘門較近位置，摻氣現(xiàn)象僅在靠近自由面一定深度內(nèi)發(fā)生，空氣還未貫穿至隧洞底部，因而底部還存在清水區(qū)。隨著與閘門距離的增加，摻氣現(xiàn)象進一步發(fā)生，摻入水體的空氣逐漸到達(dá)隧洞底部，清水區(qū)逐漸消失。此外，摻氣濃度分布曲線變化趨勢更加平緩，相比靠近閘門位置分布更加均勻。當(dāng)來流水頭逐漸增大時，各斷面處摻氣濃度分布基本一致。對于某一斷面位置而言，由于來流水頭逐漸增大，來流流速增大，水-氣紊動作用更加劇烈，摻氣量增大，造成摻氣水深增大，摻氣濃度分布更加平緩。在x=1.00 m斷面處，3種來流水頭下隧洞底部附近均存在清水區(qū)，當(dāng)來流水頭增大時，自由面位置摻氣區(qū)域深度越大，且摻氣水深越大；在x=5.98 m斷面處，隧洞底部附近清水區(qū)均已消失，隨著來流水頭增大，近底摻氣濃度增大。

圖8 閘門開度為40%時不同來流水頭He下各斷面摻氣濃度分布

當(dāng)閘門開度為40%時，不同來流水頭近底摻氣濃度沿程變化如圖9所示。由圖9可見，距離閘門較近的位置，摻入水體的空氣量較少，未擴散至隧洞底部，因而近底摻氣濃度均為零，隨著與閘門距離的增大，空氣在劇烈的紊動作用下擴散至隧洞底部，近底摻氣濃度大于0。由出口位置近底摻氣濃度可明顯得到：隨著來流水頭He的逐漸增大，水-氣紊動作用更加劇烈，造成出口位置近底摻氣濃度由He=5 m時的0.006增大至He=15 m時的0.060。

圖9 閘門開度為40%時不同來流水頭He下近底摻氣濃度沿程變化

摻氣層厚度反映了空氣對摻氣水流的影響范圍，將摻氣層厚度定義為0.01

(12)

式中，haw為摻氣層厚度，m；ZC=0.01為C=0.01時的水深，m；ZC=0.99為C=0.99時的水深，m。

閘門開度為40%時，不同來流水頭摻氣層厚度沿程變化如圖10所示。圖10中結(jié)果表明，對于不同來流水頭，摻氣層厚度變化趨勢一致，即在距閘門較近位置處的摻氣層厚度較小，隨著與閘門距離的增大，摻氣層厚度開始迅速增大，流動一段距離后，增加趨勢變緩，趨于穩(wěn)定狀態(tài)。

圖10 閘門開度為40%時不同來流水頭He下?lián)綒鈱雍穸妊爻套兓?/p>

對于某一斷面來說，隨著來流水頭的增大，水流流速增大，水-氣之間的紊動作用更加劇烈，洞內(nèi)空氣更易于摻入水體中，造成水體摻氣程度進一步加強，斷面上摻氣量隨之增大，出口位置摻氣層厚度由0.08 m增大至0.12 m。

5 結(jié) 論

數(shù)值模擬是目前研究摻氣問題的重要方法，本文基于自主開發(fā)的NEWTANK數(shù)值模型，構(gòu)建了混合-濃度數(shù)值模型用于研究水氣二相流問題，通過有限差分方法進行離散并求解，并對洞內(nèi)摻氣水流進行了數(shù)值模擬研究，主要結(jié)論如下：

(1)通過混合模型獲得混合流場和壓力場，紊流模型簡化流場內(nèi)的紊動發(fā)展過程，濃度模型用于追蹤自由液面位置并獲取摻氣水流內(nèi)部的摻氣濃度分布信息。水流流速與摻氣濃度的數(shù)值模擬結(jié)果和試驗結(jié)果吻合良好，驗證了該數(shù)值模型研究洞內(nèi)摻氣問題的可行性。

(2)距閘門較近位置，隧洞底部存在清水區(qū)，隨著與閘門距離的增加，摻氣進一步發(fā)生，空氣貫穿至底部，清水區(qū)逐漸消失，摻氣濃度分布曲線變化趨勢更加平緩。當(dāng)洞內(nèi)余幅較大時，摻氣水流的拖曳能力有限，僅能抵抗部分壓差作用帶走靠近摻氣水體的空氣，洞頂位置的空氣將在壓差作用下由隧洞出口位置進入洞內(nèi)，形成一個大的氣流循環(huán)，且洞內(nèi)上方逆向氣流速度沿水流方向由2 m/s增大至5 m/s。

(3)當(dāng)閘門開度一定，來流水頭由5 m增大至15 m時，水流流速分布相似，流速逐漸增大，水-氣間紊動作用更加劇烈，洞內(nèi)空氣更易于摻入水體中，摻氣程度進一步加強，造成靠近出口位置的近底摻氣濃度由0.006增大至0.060，摻氣層厚度由0.08 m增大至0.12 m。