基于氣體釋放模型的氣液兩相流數值模擬

， , b， ，

(武漢大學 a. 動力與機械學院， b. 水力機械過渡過程教育部重點實驗室， 湖北 武漢 430072)

兩相流動模型需要很多參數來描述，流動機理復雜，流動不穩定，實驗條件苛刻，某些參數很難測量，并且相間的結構和作用形式復雜，開展氣液兩相流動的模型研究難度很大。近年來，隨著國內學者對氣液兩相流理論研究的深入，在模型的建立及數理模型的數值求解方面取得了較大進展。

Kranenburg[1]和Martin[2]發現，用Lax-Wendroff方法能較好地捕捉激波。Zielke[3]發現，對于空化持續時間超過1.5 s的氣液兩相瞬變流，應考慮氣體釋放的因素。文獻[4-6]中在考慮氣體對波速影響的基礎上，推導了氣液兩相瞬變流的波速方程。文獻[7-8]中則對氣液兩相均質瞬變流進行了計算研究。文獻[9-10]中詳細論述了塞狀流的瞬變機理，并建立了對計算和分析含氣管路十分有效的瞬變數學模型。Martin等[11]采用漂移流模型對氣液兩相瞬變流進行計算，求解塞狀流瞬變過程的波速，并通過實驗對計算結果進行了驗證。劉竹溪等[12]在20世紀80年代初就開始了火電廠循環水系統的水錘實驗研究，提出了合理的防護措施?？紤]到流體中的自由氣體對流體瞬變過程中壓力波波速的影響，楊建東等[13]對該過程中氣體釋放規律進行了研究，并推導了相應的氣體釋放公式。

在長距離輸水管道系統中，意外斷電造成的事故停泵，或者閥門等設備出現故障而突然關閉所產生的水力過渡過程現象廣泛存在。輸水系統水錘不僅會損壞管道系統中的水力元件，而且對管道系統的安全運行造成威脅。正確分析流體瞬變過程的水力特性對輸水管道系統的設計和運行有重要的工程實踐意義，不僅可以優化系統，降低工程造價，而且保證了系統的安全運行。為了解決在長距離輸水管道計算中傳統方法不夠精確，一般氣液兩相流模型計算過程復雜并且計算量大的問題，本文中通過建立氣體釋放模型，并運用規定時間間隔的特征線法求解，推導出控制方程，最后將結果與傳統方法及實驗進行對比，在提高精度的基礎上最大限度地加快運算速度。

1 氣體釋放模型的建立

實際液體中有溶解的氣體，盡管其容積比很小，但是不應忽略。氣泡在液體中的分布及析出過程是一個極其復雜的問題，相比于單相瞬變流模型的基本假設，氣體釋放模型假定當流體中的壓力小于液體的飽和壓力時，開始發生氣體釋放。該模型假定氣體釋放集中在各計算截面上，釋放出的氣體分布在整個區域，并且氣泡體積和壓力的變化遵循等溫規律，計算時需要考慮自由氣體對壓力波波速的影響。當壓力小于汽化壓力時，蒸汽穴仍可在管道固定的計算截面上生成。

1.1 氣體釋放計算

在流體瞬變過程中，當流體壓力減至液體飽和壓力時，溶解在液體中的氣體開始釋放，并向已有的氣泡擴散。流體中氣體向氣泡擴散是一個相對緩慢的物理過程，同時在瞬變過程中只有少量的氣體釋放，所溶解的氣體濃度只有輕微的變化，但是在瞬變過程中，壓力和波傳播速度與自由氣體的孔隙度有極為密切的聯系，因此建立流體的氣體釋放數學模型對氣液兩相瞬變流的計算有重要意義。

氣體釋放速率的計算公式為

(1)

1.2 控制方程

對于氣體釋放模型，其控制方程與單相瞬變流相同，只是考慮了氣體釋放而引起的波速變化。其連續性方程和運動方程為

(2)

式中：x為空間橫坐標；H為壓力水頭；c為波速；v為水的流速；D為管道直徑；θ為傾角；f為水力損失系數。

1.3 規定時間間隔特征線法

以絕對壓頭Ha為變量，并將波速方程代入特征線轉換后的連續方程和運動方程，同時忽略流體速度的影響，則可得

(3)

(4)

將式(3)兩邊同時乘以cdt/g，并將C+的相容性方程沿著特征線積分(設網格長度為dx=cdt)，同時對水力損失項進行一階近似，則對水平管可得

式中：A為橫截面積；QP為P點流量；QR為R點流量。令

其中HaR為R點的絕對壓頭，則式(5)可表示為

(6)

式中HaP為P點的絕對壓頭。

對于特征線C+，可將其變換為

(7)

將該式沿著特征線C+積分，式(7)可整理為數值形式，即

(8)

式中：xP為P點橫坐標；xR為R點橫坐標。

(9)

定義數值計算的插值度為

(10)

同理，對于C-，可以得出其離散方程為

(11)

其中

式中：QS為S點流量；HaS為S點的絕對壓頭；xS為S點橫坐標；F(HaR,HaP)和F(HaS,HaP)分別為特征線C+和C-相容性方程中的壓力估算值。

為了求解各計算截面上下一時步的HaP、QP等，需要采用迭代算法。在聯立方程(6)、(11)計算HaP時，需要用到牛頓法。具體計算流程如圖1所示。

t—計算時長；Δt—計算的時間步長。圖1 氣體釋放模型模型計算流程圖

1.4 氣體釋放模型邊界條件

氣體釋放模型具有與單相瞬變流模型相同的邊界條件， 具體表述如下： 對于上游邊界條件， 首先利用式(10)計算插值度ζ， 然后按照標準線性插值方法求得HaS和QS， 進而求得M，而上游已知壓力變化曲線，利用特征線C-，將其代入式(11)可得第一截面的流量QS(1)，利用已求得的第一截面的壓力HP(1)和流量QP(1)重新計算插值度ζ。重復上述過程，直到得到要求的計算精度為止，即壓力為已知邊界條件，流量為計算邊界條件。

管道下游為恒定水位的水庫，已知最末端截面壓力HP(2)為常數。首先利用式(10)計算插值度ζ，然后按照標準線性插值方法求得HaR和QR，進而求得N，利用特征線C+，將N代入式(3)可得最末端截面流量QP(2)，利用已求得的HP(2)和QP(2)重新計算插值度ζ。重復上述過程，直到得到要求的計算精度為止。

2 求解實例

計算中應用了荷蘭代爾夫特理工大學水力學實驗室中水平管的上游端模擬泵的故障停車和重新開車實驗[1]。在上游端規定一個壓力-時間特性曲線來模擬泵的壓力特性，下游是一個保證壓力恒定的大水箱。水力學實驗系統裝置如圖2所示。

圖2 水力學實驗系統裝置

實驗系統參數如下：實驗管道長度L為1 450 m，直徑D為0.1 m，水力損失系數f為0.016 5，汽化壓力Hv為0.2 m，飽和蒸汽壓頭Hs為7.75 m，流量Q為0.015 8 m3/s，初始含氣率φ為5.0×10-6。氣體體積彈性模量Kg為2.068 5×105Pa，水的體積彈性模量Kl為2.18×109Pa，鋼管的楊氏模量E為2.0×1011Pa，管壁厚度e為10 mm，氣體和水的初始密度ρg、ρl分別為1.205、 998.2 kg/m3。模擬上游壓力變化曲線如圖3所示。

圖3 上游壓力-時間曲線

3 單相瞬變流模型與氣體釋放模型計算結果及分析

單相瞬變流模型和氣體釋放模型的數值方法均采用傳統特征線法，將水平管分為20個等長度的管段，用FORTRAN語言編程進行計算。其中，單相瞬變流模型結合相應的邊界條件，采用假定波速固定的特征線方程進行離散計算；氣體釋放模型結合相應的邊界條件，采用考慮自由氣體對波速影響的變波速的特征線方程(4)、(9)進行離散計算。

單相瞬變流模型計算結果與代爾夫特理工大學實驗曲線對比如圖4所示。由圖4(a)可知，單相瞬變流模型的計算結果中存在很多與實驗結果不符的壓力尖峰，計算的壓力峰值98 m比實驗中出現的壓力峰值77 m增大27.3%。圖4(b)中也存在這種現象。以上結果說明采用單相瞬變流模型進行計算會產生不真實的壓力增大。

(a)0.4 L

(b)0.8 L L—實驗管道長度1 450 m。

氣體釋放模型計算結果與實驗曲線對比如圖5所示。由圖可知，與單相瞬變流模型的計算結果相比， 壓力的變化過程更符合實驗曲線氣體釋放模型對壓力第1個峰值模擬很準確，但第2個壓力峰值的模擬結果相位相應提前，同時壓力峰值偏小。

(a)0.4 L

(b)0.8 L L—實驗管道長度1 450 m。

單相瞬變流模型計算結果與氣體釋放模型計算結果對比如圖6所示。由圖可知，氣體釋放模型計算曲線比單相瞬變流模型更光滑，2種模型在第1個波峰和波谷均能比較好地預測水錘壓力，但是第2個波峰與波谷可以較明顯地看出氣體釋放模型更能體現實際水錘的變化規律。

4 結論

本文中對荷蘭代爾夫特理工大學水力學實驗室中水平管的上游端模擬泵的故障停車實驗進行了數值模擬，比較了傳統的單相瞬變流模型和氣體釋放模型的計算結果，得出以下結論：

1)單相瞬變流模型計算的壓力峰值偏大，會產生不真實的壓力增大，原因是沒有考慮氣體會導致水錘波速減小，進而導致水錘壓力減小。

2)氣體釋放模型計算結果與實驗曲線比較符合， 能夠準確地模擬第1個壓力峰值， 但其對第2個壓力峰值的模擬結果相位相應提前， 同時壓力峰值偏小， 說明氣體釋放后， 漂移、 部分重新溶解等一系列復雜因素導致第2個壓力峰值的計算結果有所偏差。

(a)0.4 L

(b)0.8 L L—實驗管道長度1 450 m。

3)整體來看，采用氣體釋放模型能更好地反映實驗裝置中流體流動的真實特性。