三維挑流水舌及其碰撞的LBM數(shù)值模擬研究

彭燕祥，張 華，吳國華

（華北電力大學 水利與水電工程學院，北京 102206）

1 研究背景

隨著中國水電事業(yè)的不斷發(fā)展，水電站向高水頭、大流量的方向發(fā)展。泄洪霧化問題越發(fā)突出，尤其是采用挑流消能方式泄洪時，會在水電站下游局部區(qū)域形成大面積的雨霧。相關(guān)的研究表明：挑流泄洪霧化與水舌自身摻氣擴散、水舌之間相互碰撞和水舌入水噴濺密切相關(guān)，要對挑流泄洪霧化進行更深入的研究，就必須對挑流水舌的運動特性，以及水舌與空氣的耦合作用引起的水舌風場進行研究。

針對水舌的運動特性，許多學者建立不同的數(shù)學模型進行了研究，劉宣烈等［1-3］通過實驗對水舌運動過程的摻氣和擴散過程進行研究，得到了水舌運動的軌跡方程以及水舌橫、縱斷面的擴散規(guī)律。梁在潮等［4］對水舌運動軌跡、水舌摻氣量以及霧源量等參數(shù)進行了研究，并對這些參數(shù)進行擬合。張華等［5］通過建立水舌運動微分方程，并用數(shù)值方法進行求解，來研究空中水舌運動特性。劉士和等［6］運用相似理論在考慮空氣阻力及卷吸摻氣的作用下，建立了挑流水舌的數(shù)學模型。劉繼廣等［7］依據(jù)摻氣濃度相似律準則進行氣水兩相流模擬研究。張貝貝等［8］提出一種在排放口近區(qū)考慮射流垂向效應、在全計算域使用二維平面模型的計算方法。但是這些研究僅僅是對水舌進行了研究，由水舌引起的水舌風場并沒有涉及。

為了研究水舌與空氣的耦合作用引起的水舌風場。李敏等［9］采用非靜力閉合邊界層模式，模擬了水電站泄洪條件下水舌風的發(fā)展過程。彭燕祥等［10］用CFD方法模擬了二灘水電站泄洪時水舌的運動及水舌風場的空間分布。彭莘仔等［11］采用N-S方程模擬射流的水動力過程。但是這些方法都是基于有限元網(wǎng)格的歐拉方法，這個方法對于流體運動的范圍有其局限性，網(wǎng)格的劃分密度以及范圍都對誤差的影響很大。

格子玻爾茲曼方法（Lattice Boltzmann Method，LBM）是一種介觀流體系統(tǒng)模擬方法，兼具宏觀流體連續(xù)性和微觀分子動力學方法的優(yōu)點，而且計算效率較高，現(xiàn)已被廣泛應用于多組分多相流、化學反應流、氣固兩相流的研究［12-14］。在水汽兩相流兩方面，宋歌等［15］應用LBM方法對離心式霧化噴嘴出口處的液膜厚度、霧化錐角以及其他霧化特性進行了研究。葉永等［16］采用LBM對三維潰壩流進行了數(shù)值模擬與驗證。兩者都沒有對空氣流場進行研究。張華等［17］應用格子Boltzmann兩相流模型，對水舌在空氣中的運動進行了數(shù)值模擬，但是模擬的水舌是二維的，且研究的水舌尺度較小。

為了分析三維挑流水舌的運動特性，以及兩股挑流水舌的碰撞特性，本文采用LBM方法對挑流水舌進行三維數(shù)值模擬研究。

2 格子玻爾茲曼方法

粒子運動的速度矢量見圖1。

圖1 D3Q19網(wǎng)格模型與粒子速度向量

3 挑流水舌空中運動數(shù)值模擬

3.1 計算模型與參數(shù)設(shè)置 采用長江水科院［20］對泄洪霧化概化模型，進行數(shù)值模擬研究，如圖2所示，泄水建筑物長2.3 m，高1.7 m，寬0.28 m，出口挑坎挑角20°。數(shù)值模擬工況采用了典型的實驗工況，即水頭0.244 m，水頭落差2.05 m。

圖2 數(shù)值模擬幾何模型

3.2 模型參數(shù)設(shè)置 模型采用xflow兩相流計算模型，分別為液體相和空氣相，液體相密度為998.3 kg/m3，動力黏度為0.001 Pa·s?？諝庀嗝芏葹?.225 kg/m3，動力黏度為1.7894 10-5Pa·s。

由于時間步長和解析度直接關(guān)系到計算結(jié)果的精度和計算效率。所以在保證計算結(jié)果精度的同時還必須考慮計算時間。經(jīng)多方面考慮以及多次試算之后，選擇時間步長為0.001 s，解析度為0.03 m最為合適。

環(huán)境參考壓力取一個標準大氣壓，邊界條件1、2截面采用壓力進口條件，3截面采用壓力出口，相對壓力都設(shè)置為0。其余面均設(shè)置為自由滑移壁面邊界條件。

3.3 水舌風場計算結(jié)果 圖3為t=2.98 s時的水舌風場的空間分布計算圖，從圖中可以比較清楚的看出水舌與空氣運動交界面，挑流水舌隨著水舌的運動不斷的摻氣擴大。水舌的運動帶動周圍的空氣沿水舌方向運動，并且在水舌入水處，水舌表面空氣流速與水舌的流速相近風速達到6 m/s，風向沿水舌速度方向。

圖3 水舌風場的空間分布

由于水舌下緣在臨近下游水面區(qū)域與下游水面形成一個相對閉塞的空間，特別時水舌入水附近，如圖3（c）所示，水舌下緣的風場情況相較于水舌上緣相對更加的復雜，在水舌下沿出現(xiàn)較多的漩渦流場，如圖3（d）所示。從圖中可以看出水舌入水與下游水面碰撞，引起下游水面液滴飛濺，在下游水面形成一個與水舌入水面為邊界的水滴噴濺區(qū)，并且沿水舌外緣入水區(qū)域的噴濺水滴，由于水舌入水夾角的關(guān)系，沿x方向的速度分量比其他區(qū)域的大。這也是泄洪霧化的主要霧源。

3.4 水舌運動特性驗證分析 數(shù)值模擬采用了典型工況：水頭0.244 m，水頭落差2.05 m。對水舌的內(nèi)緣和外緣軌跡和挑距，以及水舌的厚度沿水舌軌跡的變化等進行分析。

采用2次多項式對水舌的內(nèi)緣、外緣以及中心內(nèi)線進行擬合，分別得到：

內(nèi)緣曲線表達式為：

外緣曲線表達式為：

中心線曲線表達式：

挑流水舌的外緣和內(nèi)緣，分別與文獻［5］做了對比，如圖4所示。

圖4 挑流水舌外緣和內(nèi)緣的空間變化

從圖4中可以看出，兩者的水舌軌跡基本一致。在水舌的中間段，本文的外緣線比文獻［5］的結(jié)果偏高，中間段水舌的擴散速度比文獻［5］更快。

圖5為水舌流速與下游水位截面位置的體積分布計算結(jié)果，水舌入水的速度在6 m/s左右。水舌內(nèi)緣挑距約為2.8 m，外緣挑距約為3.5 m，橫向?qū)挾燃s為0.75 m。與實測數(shù)據(jù)對比誤差分別為9.6%、7.9%、7.1%，如表1所示，可以看出挑距的模擬值相對實測數(shù)據(jù)偏低。水舌入水的形狀如圖5（b）所示，大致為一個圓角的矩形趨近于橢圓形，與文獻［21］中單孔泄洪實驗的水舌落水形狀基本一致，這可以為水舌入水噴濺霧源的進一步研究提供參考。

表1 水舌挑距實測值［17］與模擬值的對比

圖5 挑流水舌俯視圖與水舌入水截面圖

圖6為水舌模擬與理論計算的厚度與長度關(guān)系變化圖，其中H為水舌厚度/m，H0為水舌初始斷面厚度/m，S為水舌長度/m。通過挑流水舌厚度與水舌長度的關(guān)系圖，可以將水舌的厚度變化可以分成三個階段，如圖的Ⅰ、Ⅱ、Ⅲ階段。

圖6 挑流水舌厚度與水舌長度的關(guān)系

Ⅰ階段水舌擴散緩慢，且水舌的厚度的增長與水舌長度呈線性關(guān)系，這是由于Ⅰ階段水舌尚有核心區(qū)存在，水舌摻氣量少，水體比較密實。

Ⅱ階段水舌的厚度擴散速度先加快后慢慢變緩，這是因為由于水舌速度的增加，空氣對水舌的阻力增加，這樣水舌內(nèi)部的紊動也就隨之增加，引起水舌的破碎和大量摻氣，使水舌厚度快速的增加，隨著水舌摻氣的增加逐漸趨于飽和，使得摻氣速度逐漸趨于變緩，水舌的厚度增加速度變緩。

Ⅲ階段為水舌的均勻擴散段，水舌的厚度的增長與水舌長度呈現(xiàn)近似線性變化。這與文獻［1］中通過實驗數(shù)據(jù)分析，將水舌厚度的變化分為三個不同變化區(qū)段，即：（1）出口緩擴段；（2）水舌急劇擴散段；（3）水舌均勻擴散段；結(jié)論基本一致。

如圖7為水舌速度與水舌長度的關(guān)系對比圖，其中的文獻［5］-v、文獻［5］-vx、文獻［5］-vy為采用文獻［5］的數(shù)學模型，忽略空氣阻力情況下的計算結(jié)果。從圖中可以看出本文-v和文獻［5］-v隨著S增加先減小后增加、本文-v和文獻［5］-v的上升階段，本文-v相對與文獻［5］-v的值偏小，這是因為采用文獻［5］數(shù)學模型計算水舌時，沒有考慮空氣阻力以及水舌風場的影響。

圖7 水舌速度與水舌長度的關(guān)系

本文-vx值隨著水舌挑S的增加而減小，而文獻［5］-vx水舌由于沒有考慮空氣阻力的影響，水舌速度文獻［5］-vx的值隨著水舌挑距S的增加幾乎不變。

本文-vy和文獻［5］-vy隨S增加先減小到0后增加，當本文-vy和文獻［5］-vy等于0時，水舌上升到最高點。本文-vy與文獻［5］-vy隨S變化趨勢基本是一致，這是由于在y方向水舌雖然受重力和空氣阻力以及水舌風場的共同作用，但是重力起主導作用。

3.5 水舌渦量分析 渦量為流體速度的旋度：

渦量是描寫旋渦運動的重要物理量之一，在xy平面取水舌的中心平面得水舌的中間截面的渦量云圖，即渦量在z軸的投影用Ωz表示，如圖8所示。從圖中可以看出水舌的上下表面在初始階段Ωz遠遠高于內(nèi)部水舌，這是因為水舌初始階段摻氣并不充分，僅僅外部水舌與空氣作用，使得水舌的上下表面的Ωz遠遠高于內(nèi)部水舌。這也是水舌摻氣的重要影響因素，隨著水舌的摻氣量不斷變大，Ωz向水舌的內(nèi)部輸運，在水舌的末端Ωz達到最大，Ωz達到139 s-1。

圖8 水舌中心截面渦量云圖

4 水舌碰撞的數(shù)值模擬

格子Boltzmann方法可以更加容易的實現(xiàn)兩股水舌的碰撞，以及由于水舌碰撞而形成的相對復雜的水汽交界面，如圖9所示。

圖9 兩股水舌碰撞水舌速度分布圖

圖9和圖10分別為水舌碰撞的數(shù)值模擬圖與文獻［22］中水舌碰撞的實驗圖象。從圖中可以看出兩者的碰撞水舌形態(tài)基本一致，水舌經(jīng)過碰撞后，水花四濺，水舌在碰撞后形成一個錐形的射流形態(tài)，并且由于兩股水舌的碰撞形成的壓力波影響，使得上股水舌的下表面破碎形成一股反向微射流，微射流以液滴為主且運動狀態(tài)較為復雜，但其大多數(shù)液滴的運動方向與水舌的運動方向相反。這為進一步研究多股水舌泄流，以及水舌入水噴濺等復雜情況提供參考。

圖10 兩股水舌碰撞的試驗結(jié)果圖［22］

5 結(jié)論

本文采用格子Boltzmann方法針對三維挑流水舌以及兩股水舌的碰撞進行數(shù)值模擬研究，得到如下結(jié)論：

（1）格子Boltzmann方法能很好的模擬泄洪挑流水舌和水舌風場，水舌的運動界面清晰可見。并且不用對水舌的周圍的大片區(qū)域進行網(wǎng)格劃分，計算更加便捷，效率更高。水舌內(nèi)緣挑距約為2.8 m，外緣挑距約為3.5 m，橫向?qū)挾燃s為0.75 m。與實測數(shù)據(jù)對比誤差分別為9.6%、7.9%、7.1%，挑距的模擬值相對實測數(shù)據(jù)偏低。為泄洪挑流水舌的流場研究提供新的方法。

（2）挑流水舌在初始階段由于水舌與空氣的相互作用，使得水舌的上下表面的Ωz遠遠高于內(nèi)部水舌，這也是水舌摻氣的重要影響因素，隨著水舌的摻氣量不斷變大，Ωz向水舌的內(nèi)部輸運，在水舌的末端渦量達到最大，渦量達到139 s-1。

（3）格子Boltzmann方法可以對兩股水舌的碰撞進行模擬，水舌碰撞的數(shù)值模擬與實驗條件下水舌運動形態(tài)基本一致，水舌經(jīng)過碰撞后，水花四濺，水舌在碰撞后形成一個錐形的射流形態(tài)，并且由于兩股水舌的碰撞形成的壓力波影響，使得上股水舌的下表面破碎形成一股反向微射流。為多孔泄流的復雜情況，水舌的運動特性以及流場研究提供參考。