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

， , b， ，

(武漢大學(xué) a. 動(dòng)力與機(jī)械學(xué)院， b. 水力機(jī)械過渡過程教育部重點(diǎn)實(shí)驗(yàn)室， 湖北 武漢 430072)

兩相流動(dòng)模型需要很多參數(shù)來描述，流動(dòng)機(jī)理復(fù)雜，流動(dòng)不穩(wěn)定，實(shí)驗(yàn)條件苛刻，某些參數(shù)很難測量，并且相間的結(jié)構(gòu)和作用形式復(fù)雜，開展氣液兩相流動(dòng)的模型研究難度很大。近年來，隨著國內(nèi)學(xué)者對氣液兩相流理論研究的深入，在模型的建立及數(shù)理模型的數(shù)值求解方面取得了較大進(jìn)展。

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

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

1 氣體釋放模型的建立

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

1.1 氣體釋放計(jì)算

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

氣體釋放速率的計(jì)算公式為

(1)

1.2 控制方程

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

(2)

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

1.3 規(guī)定時(shí)間間隔特征線法

以絕對壓頭Ha為變量，并將波速方程代入特征線轉(zhuǎn)換后的連續(xù)方程和運(yùn)動(dòng)方程，同時(shí)忽略流體速度的影響，則可得

(3)

(4)

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

式中：A為橫截面積；QP為P點(diǎn)流量；QR為R點(diǎn)流量。令

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

(6)

式中HaP為P點(diǎn)的絕對壓頭。

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

(7)

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

(8)

式中：xP為P點(diǎn)橫坐標(biāo)；xR為R點(diǎn)橫坐標(biāo)。

(9)

定義數(shù)值計(jì)算的插值度為

(10)

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

(11)

其中

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

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

t—計(jì)算時(shí)長；Δt—計(jì)算的時(shí)間步長。圖1 氣體釋放模型模型計(jì)算流程圖

1.4 氣體釋放模型邊界條件

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

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

2 求解實(shí)例

計(jì)算中應(yīng)用了荷蘭代爾夫特理工大學(xué)水力學(xué)實(shí)驗(yàn)室中水平管的上游端模擬泵的故障停車和重新開車實(shí)驗(yàn)[1]。在上游端規(guī)定一個(gè)壓力-時(shí)間特性曲線來模擬泵的壓力特性，下游是一個(gè)保證壓力恒定的大水箱。水力學(xué)實(shí)驗(yàn)系統(tǒng)裝置如圖2所示。

圖2 水力學(xué)實(shí)驗(yàn)系統(tǒng)裝置

實(shí)驗(yàn)系統(tǒng)參數(shù)如下：實(shí)驗(yàn)管道長度L為1 450 m，直徑D為0.1 m，水力損失系數(shù)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 上游壓力-時(shí)間曲線

3 單相瞬變流模型與氣體釋放模型計(jì)算結(jié)果及分析

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

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

(a)0.4 L

(b)0.8 L L—實(shí)驗(yàn)管道長度1 450 m。

氣體釋放模型計(jì)算結(jié)果與實(shí)驗(yàn)曲線對比如圖5所示。由圖可知，與單相瞬變流模型的計(jì)算結(jié)果相比， 壓力的變化過程更符合實(shí)驗(yàn)曲線氣體釋放模型對壓力第1個(gè)峰值模擬很準(zhǔn)確，但第2個(gè)壓力峰值的模擬結(jié)果相位相應(yīng)提前，同時(shí)壓力峰值偏小。

(a)0.4 L

(b)0.8 L L—實(shí)驗(yàn)管道長度1 450 m。

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

4 結(jié)論

本文中對荷蘭代爾夫特理工大學(xué)水力學(xué)實(shí)驗(yàn)室中水平管的上游端模擬泵的故障停車實(shí)驗(yàn)進(jìn)行了數(shù)值模擬，比較了傳統(tǒng)的單相瞬變流模型和氣體釋放模型的計(jì)算結(jié)果，得出以下結(jié)論：

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

2)氣體釋放模型計(jì)算結(jié)果與實(shí)驗(yàn)曲線比較符合， 能夠準(zhǔn)確地模擬第1個(gè)壓力峰值， 但其對第2個(gè)壓力峰值的模擬結(jié)果相位相應(yīng)提前， 同時(shí)壓力峰值偏小， 說明氣體釋放后， 漂移、 部分重新溶解等一系列復(fù)雜因素導(dǎo)致第2個(gè)壓力峰值的計(jì)算結(jié)果有所偏差。

(a)0.4 L

(b)0.8 L L—實(shí)驗(yàn)管道長度1 450 m。

3)整體來看，采用氣體釋放模型能更好地反映實(shí)驗(yàn)裝置中流體流動(dòng)的真實(shí)特性。