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

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

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

0 引 言

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

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

本文利用基于格子玻爾茲曼方法的新型無網格CFD方法分析三維重力流管道系統的充水和排空過程，以新的數值計算方法分析系統級別的流場瞬變過程是一個創新的研究方向。為節約計算資源和時間，僅取長度較短的小型管道系統作為研究對象，著重從原理上說明計算的可行性。

1 數學模型

1.1 重力流數學模型

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

圖1 常見重力流系統Fig.1 Common system of gravity flow

在重力流流動過程中，整個系統機械能守恒，有：

(1)

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

當管流為滿管流時，由于管徑和密度不變，則進、出口質量流量相同，流速相同，即：

m1=m2=m

(2)

u1=u2=u

(3)

此時，系統的能量損失包括入口損失、管道損失和出口損失，E的表達式為：

(4)

式中：k1和k2分別是管道進口和出口的阻力系數；D是管道直徑，m；Li是不包括進出口的分段管道當量長度，m；λi是分段管道沿阻系數；ζj是局阻系數。

以水頭形式表示重力流的總推動力為：

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

(5)

1.2 格子玻爾茲曼方法

在物理學中對流體流動描述有三個層次：宏觀、微觀和介觀[16]?；诮橛^模型的格子玻爾茲曼方法以流體的分子運動論為基礎，根據微觀運動過程的基本特征建立了簡化的時間和空間完全離散的動力學格子模型[17]：在格子之間分布的眾多介觀流體粒子根據一定的法則沿格線遷移并在格點上發生碰撞，運動過程遵從動力學定律并服從統計定律，通過對流場中大量粒子的統計平均得出流體的宏觀運動變量。

格子玻爾茲曼方法與傳統CFD的N-S方程算法相比，最大的特點就是求解的無網格化，只需定義格子尺寸即可。另外其還有許多獨特的優點：①從微觀粒子角度進行模擬，可以比較方便地處理不同組分、不同流體間界面間的復雜相互作用不需要借助經驗和半經驗公式；②相對于N-S方程的非線性對流項，格子玻爾茲曼方法相空間對流過程是線性的，求解更容易，演化過程更清晰；③壓力用狀態方程表示，不必與速度迭代求解，算法比較簡單。

采用單松弛時間算子[18]且帶體積力項的格子玻爾茲曼方程如下：

(6)

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

碰撞步：

(7)

遷移步：

(8)

而局部平衡態分布函數為：

(9)

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

2 計算方案及結果

2.1 重力流管道充水

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

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

重力流管道充水過程的管內液面變化如圖3所示，其中管道充滿水的時間為95 s。

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

管道進口閘閥開啟過程的管內流速變化如圖4所示。

圖4 開閥充水過程管內速度場變化情況Fig.4 The changing velocity distribution with the opening gate valve in filling process

取位于進口閘閥后面1m處管軸線上的一點作為監測點，提取重力流管道充水過程中該點的壓力和流速值，變化曲線分別如圖5和6所示。

圖5 充水過程監測點壓力變化Fig.5 The changing pressure of monitoring point in filling process

圖6 充水過程監測點流速變化Fig.6 The changing velocity of monitoring point in filling process

監測點位于閥后軸線上，受開閥產生的紊流影響，該點處的流速和壓力有較大波動。閥門開度較小時，水流主要沿管壁流動，監測點上的壓力和流速變化較小。而閘閥開啟50s左右時，壓力和流速均有一快速上升過程，隨后壓力持續下降，而流速則維持比較穩定的狀態。

2.2 重力流管道排空

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

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

重力流管道排空過程的管內液面變化如圖8所示，管道完全排空的時間為89 s。

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

管道出口閘閥開啟過程的流速變化如圖9所示。

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

圖9 開閥排空過程管內速度場變化情況Fig.9 The changing velocity distribution with the opening gate valve in draining process

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

圖10 排空過程監測點壓力變化Fig.10 The changing pressure of monitoring point in draining process

圖11 排空過程監測點流速變化Fig.11 The changing velocity of monitoring point in draining process

出口閘閥離管道出口較近，受閥門啟動擾動，管道出口監測點在排空過程中產生了一定的負壓，但在閥門完全打開(60 s)一段時間后，在80 s左右出現快速上升并回歸正值。另外，整個排空過程中管道出口流速小幅波動地上升，在89 s時流速驟降至零，說明此時管內液體已完全排空。

2.3 重力流充排水一維計算結果

為驗證重力流充水和排空過程的三維瞬態計算結果，對該管道系統同時進行一維系統級別的對比計算?；谔卣骶€法的一維管道計算不能得到管內流場變化情況，而著重求得管道重力流充水和排空過程閥后壓力、流速等參數的整體變化情況，并與三維計算結果進行對比驗證。管道重力流充排水一維計算模型如圖12所示。

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

提取與三維計算對應的充水和排空過程閥后監測點的壓力和流速值，變化曲線如圖13和圖14所示。

圖13 充水過程監測點壓力、流速變化(1D)Fig.13 The changing pressure and velocity of monitoring point in filling process(1D)

圖14 排空過程監測點壓力、流速變化(1D)Fig.14 The changing pressure and velocity of monitoring point in draining process(1D)

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

在滿管排空階段，三維與一維仿真相比，監測點壓力波動較大；而兩者流量變化趨勢相同——上升至峰值后在89 s左右迅速下降至零，表明兩個模型的排空時間相同。

另外，在管道重力流充水和排空的一維仿真中，壓力和流速變化曲線的連續性較好，因其只考慮管軸線上的參數變化；而對應的三維瞬態模擬中參數變化過程則波動較大，體現了三維空間變化的影響。

3 總 結

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

(1)管道充水時管內液面變化過程和充水時間(95 s)；

(2)管道充水時進口閘閥開啟過程的流態變化；

(3)管道充水時閥后監測點的壓力和流速變化曲線；

(4)管道排空時管內液面變化過程和排空時間(89 s)；

(5)管道排空時出口閘閥開啟過程的流態變化；

(6)管道排空時管道出口監測點的壓力和流速變化曲線。

同時對相同的管道系統進行基于特征線法的一維計算，并與三維模擬結果進行對比和驗證。將傳統的管道系統研究由一維分析轉變為三維的CFD分析是一種有益嘗試，結合新型的格子玻爾茲曼算法，使流場仿真結果精度更高。

□

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

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

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

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

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

[6] 甘延標. 多相流系統的格子玻爾茲曼方法與模擬研究[D]. 北京:中國礦業大學(北京), 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方法的建筑流場仿真研究[D]. 西安:西安建筑科技大學, 2010.

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

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

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

[12] 許友生, 李華兵, 方海平, 等. 用格子玻爾茲曼方法研究流動-反應耦合的非線性滲流問題[J]. 物理學報, 2004,(3):773-777.

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

[14] 李 強, 余 凱, 宣益民. 納米流體流動與傳熱過程的LBM并行計算[J]. 南京理工大學學報(自然科學版), 2005,29(6):631-634.

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

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

[17] 段欣悅. 格子玻爾茲曼方法的理論研究與應用[D]. 青島: 中國石油大學(華東), 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] 陶 實. 格子玻爾茲曼方法在計算流體力學中的應用研究[D]. 遼寧大連:大連理工大學, 2013.