基于格子玻爾茲曼的管道重力流充水和排空三維瞬變過(guò)程模擬

孟兆芳，蔣 勁，李燕輝，符向前

(1. 甘肅省水利廳水利管理局, 蘭州 730000; 2. 武漢大學(xué)動(dòng)力與機(jī)械學(xué)院, 武漢 430072)

0 引 言

常見(jiàn)的壓力管道輸水有泵流和重力流兩種輸送方式。泵流輸水系統(tǒng)主要包括泵站、輸水管道和其他水工建筑物等，以水泵作為系統(tǒng)能量源，將用水輸送到遠(yuǎn)方。而重力流系統(tǒng)則完全依賴水流的落差重力作用，不需要水泵加壓，將水從水源地輸送到遠(yuǎn)方用水戶。重力流輸水方式最大限度利用了水流的落差勢(shì)能，運(yùn)行成本低，經(jīng)濟(jì)性好，是理想的輸水方式[1]。隨著經(jīng)濟(jì)建設(shè)的快速發(fā)展，工業(yè)和人們生活用水量迅速提高，各地大中型重力流管道輸水工程也越來(lái)越多。保證重力流管道安全運(yùn)行是一道工程難題，包括流量控制、水錘壓力防護(hù)和進(jìn)排氣等眾多問(wèn)題[2]。重力流的常見(jiàn)研究方法是利用一維系統(tǒng)分析軟件進(jìn)行數(shù)值模擬，但當(dāng)管徑較大時(shí)，徑向流場(chǎng)變化明顯，將管道空間簡(jiǎn)化為一維的管軸線會(huì)使模擬結(jié)果失真。另外，一維仿真對(duì)系統(tǒng)內(nèi)的泵閥啟閉過(guò)程和管道轉(zhuǎn)折處等特殊部位的瞬態(tài)流場(chǎng)分析也不適用。

格子玻爾茲曼方法是一種新興的流體數(shù)值計(jì)算方法，其形式簡(jiǎn)單，計(jì)算局部性強(qiáng)，能靈活處理復(fù)雜的邊界問(wèn)題，且完全并行，非常適合在大規(guī)模計(jì)算機(jī)群上運(yùn)行[3]。目前該方法已在國(guó)內(nèi)外被廣泛應(yīng)用于流動(dòng)和傳熱領(lǐng)域的研究，包括：對(duì)湍流、多相流等基礎(chǔ)流體力學(xué)問(wèn)題進(jìn)行深入研究[4-6]；對(duì)室內(nèi)通風(fēng)換熱、房屋外風(fēng)壓和地下巷道熱流分布等建筑物流場(chǎng)進(jìn)行分析[7-9]；對(duì)彎道水流、鼓泡床離散顆粒和間隙滲流等具體工程問(wèn)題進(jìn)行模擬[10-12]；對(duì)流動(dòng)與傳熱過(guò)程進(jìn)行微尺度的探討等[13,14]。

本文利用基于格子玻爾茲曼方法的新型無(wú)網(wǎng)格CFD方法分析三維重力流管道系統(tǒng)的充水和排空過(guò)程，以新的數(shù)值計(jì)算方法分析系統(tǒng)級(jí)別的流場(chǎng)瞬變過(guò)程是一個(gè)創(chuàng)新的研究方向。為節(jié)約計(jì)算資源和時(shí)間，僅取長(zhǎng)度較短的小型管道系統(tǒng)作為研究對(duì)象，著重從原理上說(shuō)明計(jì)算的可行性。

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

1.1 重力流數(shù)學(xué)模型

圖1為常見(jiàn)的重力流系統(tǒng)[15]，液體從V-1截面流至V-2截面。其中，V-1處的壓力和溫度分別為P0(Pa)和T0(K)，管道入口處為P1和T1，V-2(管道出口水平面)處為P2和T2。以地面為基準(zhǔn)，三處高度分別是H0(m)、H1和H2。液體在高差的驅(qū)動(dòng)下由V-1流向V-2。為方便分析，對(duì)系統(tǒng)做出一些假定：管徑無(wú)變化，液體密度恒定，流動(dòng)過(guò)程無(wú)閃蒸等復(fù)雜現(xiàn)象，同時(shí)系統(tǒng)不涉及傳熱過(guò)程，即T0、T1和T2相等。

圖1 常見(jiàn)重力流系統(tǒng)Fig.1 Common system of gravity flow

在重力流流動(dòng)過(guò)程中，整個(gè)系統(tǒng)機(jī)械能守恒，有：

(1)

式中：m1和m2分別是進(jìn)口和出口質(zhì)量流量，kg/s；u1和u2分別是進(jìn)口和出口流速，m/s；ρ是液體密度，kg/m3；E為系統(tǒng)能量損失，J。

當(dāng)管流為滿管流時(shí)，由于管徑和密度不變，則進(jìn)、出口質(zhì)量流量相同，流速相同，即：

m1=m2=m

(2)

u1=u2=u

(3)

此時(shí)，系統(tǒng)的能量損失包括入口損失、管道損失和出口損失，E的表達(dá)式為：

(4)

式中：k1和k2分別是管道進(jìn)口和出口的阻力系數(shù)；D是管道直徑，m；Li是不包括進(jìn)出口的分段管道當(dāng)量長(zhǎng)度，m；λi是分段管道沿阻系數(shù)；ζj是局阻系數(shù)。

以水頭形式表示重力流的總推動(dòng)力為：

X=(H1-H2)+(P1-P2)/ρg

(5)

1.2 格子玻爾茲曼方法

在物理學(xué)中對(duì)流體流動(dòng)描述有三個(gè)層次：宏觀、微觀和介觀[16]?；诮橛^模型的格子玻爾茲曼方法以流體的分子運(yùn)動(dòng)論為基礎(chǔ)，根據(jù)微觀運(yùn)動(dòng)過(guò)程的基本特征建立了簡(jiǎn)化的時(shí)間和空間完全離散的動(dòng)力學(xué)格子模型[17]：在格子之間分布的眾多介觀流體粒子根據(jù)一定的法則沿格線遷移并在格點(diǎn)上發(fā)生碰撞，運(yùn)動(dòng)過(guò)程遵從動(dòng)力學(xué)定律并服從統(tǒng)計(jì)定律，通過(guò)對(duì)流場(chǎng)中大量粒子的統(tǒng)計(jì)平均得出流體的宏觀運(yùn)動(dòng)變量。

格子玻爾茲曼方法與傳統(tǒng)CFD的N-S方程算法相比，最大的特點(diǎn)就是求解的無(wú)網(wǎng)格化，只需定義格子尺寸即可。另外其還有許多獨(dú)特的優(yōu)點(diǎn)：①?gòu)奈⒂^粒子角度進(jìn)行模擬，可以比較方便地處理不同組分、不同流體間界面間的復(fù)雜相互作用不需要借助經(jīng)驗(yàn)和半經(jīng)驗(yàn)公式；②相對(duì)于N-S方程的非線性對(duì)流項(xiàng)，格子玻爾茲曼方法相空間對(duì)流過(guò)程是線性的，求解更容易，演化過(guò)程更清晰；③壓力用狀態(tài)方程表示，不必與速度迭代求解，算法比較簡(jiǎn)單。

采用單松弛時(shí)間算子[18]且?guī)w積力項(xiàng)的格子玻爾茲曼方程如下：

(6)

上式通常分解為碰撞步和遷移步兩種演化形式[19]：

碰撞步：

(7)

遷移步：

(8)

而局部平衡態(tài)分布函數(shù)為：

(9)

作為計(jì)算流體力學(xué)中一種高效的介觀數(shù)值方法，利用格子玻爾茲曼方法模擬三維重力流能得到較高的計(jì)算精度。

2 計(jì)算方案及結(jié)果

2.1 重力流管道充水

圖2是管道重力流充水計(jì)算模型，其中高位水池液位高差是50 m，管道高差60 m，管徑1 m。在管道進(jìn)口處設(shè)置一閘閥，充水時(shí)設(shè)置閘閥在60 s內(nèi)線性開(kāi)啟，總仿真時(shí)間為110 s。

圖2 重力流管道充水計(jì)算模型Fig.2 Computing model of filling gravity flowed pipes

重力流管道充水過(guò)程的管內(nèi)液面變化如圖3所示，其中管道充滿水的時(shí)間為95 s。

圖3 充水過(guò)程液面位置變化情況Fig.3 The changing liquid level in filling process

管道進(jìn)口閘閥開(kāi)啟過(guò)程的管內(nèi)流速變化如圖4所示。

圖4 開(kāi)閥充水過(guò)程管內(nèi)速度場(chǎng)變化情況Fig.4 The changing velocity distribution with the opening gate valve in filling process

取位于進(jìn)口閘閥后面1m處管軸線上的一點(diǎn)作為監(jiān)測(cè)點(diǎn)，提取重力流管道充水過(guò)程中該點(diǎn)的壓力和流速值，變化曲線分別如圖5和6所示。

圖5 充水過(guò)程監(jiān)測(cè)點(diǎn)壓力變化Fig.5 The changing pressure of monitoring point in filling process

圖6 充水過(guò)程監(jiān)測(cè)點(diǎn)流速變化Fig.6 The changing velocity of monitoring point in filling process

監(jiān)測(cè)點(diǎn)位于閥后軸線上，受開(kāi)閥產(chǎn)生的紊流影響，該點(diǎn)處的流速和壓力有較大波動(dòng)。閥門開(kāi)度較小時(shí)，水流主要沿管壁流動(dòng)，監(jiān)測(cè)點(diǎn)上的壓力和流速變化較小。而閘閥開(kāi)啟50s左右時(shí)，壓力和流速均有一快速上升過(guò)程，隨后壓力持續(xù)下降，而流速則維持比較穩(wěn)定的狀態(tài)。

2.2 重力流管道排空

圖7是重力流管道排空計(jì)算模型，初始時(shí)刻管道充滿水，液位高差60 m，管徑1 m。在管道出口處設(shè)置一閘閥，排空時(shí)設(shè)置閘閥在60 s內(nèi)線性開(kāi)啟，總仿真時(shí)間為100 s。

圖7 重力流管道排空計(jì)算模型Fig.7 Computing model of draining gravity flowed pipes

重力流管道排空過(guò)程的管內(nèi)液面變化如圖8所示，管道完全排空的時(shí)間為89 s。

圖8 排空過(guò)程液面位置變化情況Fig.8 The changing liquid level in draining process

管道出口閘閥開(kāi)啟過(guò)程的流速變化如圖9所示。

在將超高速動(dòng)能武器打擊效應(yīng)等效成淺埋爆炸后，可根據(jù)巖石中淺埋爆炸效應(yīng)，計(jì)算出超高速動(dòng)能武器打擊時(shí)，地沖擊應(yīng)力波形參數(shù)。目前計(jì)算巖石中爆炸的應(yīng)力波參數(shù)常用的計(jì)算公式為[17-18]

圖9 開(kāi)閥排空過(guò)程管內(nèi)速度場(chǎng)變化情況Fig.9 The changing velocity distribution with the opening gate valve in draining process

取位管道出口處管軸線上的一點(diǎn)作為監(jiān)測(cè)點(diǎn)，提取重力流管道排空過(guò)程中該點(diǎn)的壓力和流速值，變化曲線分別如圖10和11所示。

圖10 排空過(guò)程監(jiān)測(cè)點(diǎn)壓力變化Fig.10 The changing pressure of monitoring point in draining process

圖11 排空過(guò)程監(jiān)測(cè)點(diǎn)流速變化Fig.11 The changing velocity of monitoring point in draining process

出口閘閥離管道出口較近，受閥門啟動(dòng)擾動(dòng)，管道出口監(jiān)測(cè)點(diǎn)在排空過(guò)程中產(chǎn)生了一定的負(fù)壓，但在閥門完全打開(kāi)(60 s)一段時(shí)間后，在80 s左右出現(xiàn)快速上升并回歸正值。另外，整個(gè)排空過(guò)程中管道出口流速小幅波動(dòng)地上升，在89 s時(shí)流速驟降至零，說(shuō)明此時(shí)管內(nèi)液體已完全排空。

2.3 重力流充排水一維計(jì)算結(jié)果

為驗(yàn)證重力流充水和排空過(guò)程的三維瞬態(tài)計(jì)算結(jié)果，對(duì)該管道系統(tǒng)同時(shí)進(jìn)行一維系統(tǒng)級(jí)別的對(duì)比計(jì)算?；谔卣骶€法的一維管道計(jì)算不能得到管內(nèi)流場(chǎng)變化情況，而著重求得管道重力流充水和排空過(guò)程閥后壓力、流速等參數(shù)的整體變化情況，并與三維計(jì)算結(jié)果進(jìn)行對(duì)比驗(yàn)證。管道重力流充排水一維計(jì)算模型如圖12所示。

圖12 重力流管道一維計(jì)算模型Fig.12 1D computing model of gravity flowed pipes

提取與三維計(jì)算對(duì)應(yīng)的充水和排空過(guò)程閥后監(jiān)測(cè)點(diǎn)的壓力和流速值，變化曲線如圖13和圖14所示。

圖13 充水過(guò)程監(jiān)測(cè)點(diǎn)壓力、流速變化(1D)Fig.13 The changing pressure and velocity of monitoring point in filling process(1D)

圖14 排空過(guò)程監(jiān)測(cè)點(diǎn)壓力、流速變化(1D)Fig.14 The changing pressure and velocity of monitoring point in draining process(1D)

在空管充水階段，三維計(jì)算的監(jiān)測(cè)點(diǎn)壓力在開(kāi)閥50 s后快速下降，而一維計(jì)算的監(jiān)測(cè)點(diǎn)壓力在上升后保持在45 m水頭左右，表明三維仿真中高位水池液面下降較一維計(jì)算快。另外，兩種計(jì)算模型的監(jiān)測(cè)點(diǎn)流速變化趨勢(shì)相同，均是“快速上升-緩慢下降-平穩(wěn)”的過(guò)程，兩者流速最后均保持在16 m/s左右。

在滿管排空階段，三維與一維仿真相比，監(jiān)測(cè)點(diǎn)壓力波動(dòng)較大；而兩者流量變化趨勢(shì)相同——上升至峰值后在89 s左右迅速下降至零，表明兩個(gè)模型的排空時(shí)間相同。

另外，在管道重力流充水和排空的一維仿真中，壓力和流速變化曲線的連續(xù)性較好，因其只考慮管軸線上的參數(shù)變化；而對(duì)應(yīng)的三維瞬態(tài)模擬中參數(shù)變化過(guò)程則波動(dòng)較大，體現(xiàn)了三維空間變化的影響。

3 總 結(jié)

本文通過(guò)無(wú)網(wǎng)格的格子玻爾茲曼CFD方法模擬三維管道重力流的充水和排空過(guò)程，分別得出：

(1)管道充水時(shí)管內(nèi)液面變化過(guò)程和充水時(shí)間(95 s)；

(2)管道充水時(shí)進(jìn)口閘閥開(kāi)啟過(guò)程的流態(tài)變化；

(3)管道充水時(shí)閥后監(jiān)測(cè)點(diǎn)的壓力和流速變化曲線；

(4)管道排空時(shí)管內(nèi)液面變化過(guò)程和排空時(shí)間(89 s)；

(5)管道排空時(shí)出口閘閥開(kāi)啟過(guò)程的流態(tài)變化；

(6)管道排空時(shí)管道出口監(jiān)測(cè)點(diǎn)的壓力和流速變化曲線。

同時(shí)對(duì)相同的管道系統(tǒng)進(jìn)行基于特征線法的一維計(jì)算，并與三維模擬結(jié)果進(jìn)行對(duì)比和驗(yàn)證。將傳統(tǒng)的管道系統(tǒng)研究由一維分析轉(zhuǎn)變?yōu)槿S的CFD分析是一種有益嘗試，結(jié)合新型的格子玻爾茲曼算法，使流場(chǎng)仿真結(jié)果精度更高。

□

[1] 李 輝. 重力流壓力管道輸水系統(tǒng)的水力學(xué)問(wèn)題[J]. 山西建筑，2004,(8):73-74.

[2] 楊照明，張金承. 長(zhǎng)距離大落差重力流管道輸水技術(shù)[J]. 水利建設(shè)與管理，2001,(4):12-15.

[3] 趙 鵬, 張丹丹, 汪魯兵, 等. 格子Boltzmann并行程序的優(yōu)化與性能分析[J]. 微電子學(xué)與計(jì)算機(jī), 2008,25(10):185-188.

[4] 李 劍, 施 娟, 邱 冰, 等. 用格子玻爾茲曼方法研究粒子在渦流中的運(yùn)動(dòng)[J]. 桂林電子科技大學(xué)學(xué)報(bào), 2007,(5):387-390.

[5] 劉祖斌, 趙 鵬. 結(jié)合大渦模擬的格子玻爾茲曼方法模擬高雷諾數(shù)流動(dòng)[J]. 船舶力學(xué), 2015,(5):484-492.

[6] 甘延標(biāo). 多相流系統(tǒng)的格子玻爾茲曼方法與模擬研究[D]. 北京:中國(guó)礦業(yè)大學(xué)(北京), 2012.

[7] Crouse B, Krafczyk M, Kuhner S, et al. Indoor air flow analysis based on lattice Boltzmann methods[J]. Energy and Buildings, 2002,(34):941-949.

[8] 王 輝. 基于格子Boltzmann方法的建筑流場(chǎng)仿真研究[D]. 西安:西安建筑科技大學(xué), 2010.

[9] 韋獻(xiàn)剛. 基于格子Boltzmann方法的地下巷道熱流態(tài)仿真研究[D]. 西安:西安建筑科技大學(xué), 2009.

[10] 丁全林, 王玲玲, 汪德?tīng)? 等. 基于多松弛格子玻爾茲曼的彎道水流三維數(shù)值模擬[J]. 水科學(xué)進(jìn)展, 2012,(4):523-528.

[11] 張 博, 王利民, 王小偉, 等. 基于格子玻爾茲曼方法的單孔射流鼓泡床的離散顆粒模擬[J]. 科學(xué)通報(bào), 2013,(2):158-169.

[12] 許友生, 李華兵, 方海平, 等. 用格子玻爾茲曼方法研究流動(dòng)-反應(yīng)耦合的非線性滲流問(wèn)題[J]. 物理學(xué)報(bào), 2004,(3):773-777.

[13] 馮蕾蕾. 基于格子Boltzmann方法的復(fù)合材料微尺度傳熱特征研究[D]. 蘭州: 蘭州大學(xué), 2009.

[14] 李 強(qiáng), 余 凱, 宣益民. 納米流體流動(dòng)與傳熱過(guò)程的LBM并行計(jì)算[J]. 南京理工大學(xué)學(xué)報(bào)(自然科學(xué)版), 2005,29(6):631-634.

[15] 劉新偉. 重力流管道管徑設(shè)計(jì)[J]. 化工設(shè)計(jì), 2004,(2):20-22.

[16] 陶文銓. 數(shù)值傳熱學(xué)[M]. 2版. 西安：西安交通大學(xué)出版社, 2001.

[17] 段欣悅. 格子玻爾茲曼方法的理論研究與應(yīng)用[D]. 青島: 中國(guó)石油大學(xué)(華東), 2006.

[18] Qian Y H, Humieres D D, Lalleman P. Lattice BGK models for Navier-Stokes equation[J]. Europhysisc Letters, 1992,17(6):479-484.

[19] 陶 實(shí). 格子玻爾茲曼方法在計(jì)算流體力學(xué)中的應(yīng)用研究[D]. 遼寧大連:大連理工大學(xué), 2013.