黃 娜，周云龍，高 聚

(1.華北電力大學(xué)能源動力與機械工程學(xué)院，北京102206;2.東北電力大學(xué)能源與動力工程學(xué)院，吉林吉林132012;3.中機國能電力工程有限公司，上海200061)

0 引言

氣液兩相流換熱系統(tǒng)由于具備高效的傳熱能力，較小的溫度差分布和對抽運功率要求低等優(yōu)點，已經(jīng)在航空航天和高海拔等實際應(yīng)用中起到越來越重要的作用［1］。美國宇航局預(yù)計在未來啟動兩相流換熱系統(tǒng)來擔(dān)任主動熱控制、生命支持、推進劑輸送和原位資源利用等航空任務(wù)［2］。另外，新型兩相流循環(huán)熱控制系統(tǒng)將負責(zé)未來載人探月和探索火星工程中部分重力環(huán)境(partial gravity，10－2g0＜g ＜g0)下的能量輸運工作［3］。同時，部分重力下的氣液兩相流動也是月球基地中許多問題都涉及到的關(guān)鍵技術(shù)問題［4］。由于部分重力環(huán)境下氣液兩相的流動特性對氣液兩相流設(shè)備的工作效率具有重要影響，因此，部分重力環(huán)境下的氣液兩相流動問題成為實際應(yīng)用面臨的一個重要挑戰(zhàn)。

近年來，人們利用落管或落塔、駕駛飛機沿拋物線軌跡飛行、高空氣球搭載、火箭搭載、空間軌道飛行器搭載以及地面模擬等方法制造失重環(huán)境，進行了一些失重環(huán)境下的氣液兩相流動研究。1999年8 月，中國科學(xué)院國家微重力實驗室與俄羅斯Keldysh 研究中心展開合作，以“和平號”空間站為平臺，首次在長期、穩(wěn)定的微重力環(huán)境下進行不同重力環(huán)境下的氣液兩相流實驗研究［5］。然而，由于受到空間實驗時間的限制，在部分重力環(huán)境下僅進行了14 組實驗，難以獲取大量實驗數(shù)據(jù)來建立確定性結(jié)論。

由于長時間、穩(wěn)定的部分重力環(huán)境難以實現(xiàn)，且需要付出昂貴的代價，越來越多的學(xué)者們采用理論分析和數(shù)值模擬方法進行部分重力環(huán)境下的研究工作。趙建福等［6］以“和平號”空間站為背景研究了部分重力環(huán)境(0.1g0和0.014g0)下的氣液兩相流動特征，重點關(guān)注于流型和流型轉(zhuǎn)變界限。Clarke 和Rezkallah 等［7］基于VOF 方法建立了二維非穩(wěn)態(tài)模型來分析微重力環(huán)境下豎直管內(nèi)泡狀流的漂移速度，研究表明，液體Reynolds 數(shù)、管徑及氣泡的直徑和徑向位置對漂移速度均有重要影響。為了跟蹤豎直管內(nèi)的氣彈運動，Gabriel［8］研發(fā)了一個二維軸對稱彈狀流模型，并對μg0和g0環(huán)境下彈狀流的流場分布進行了數(shù)值研究。除此之外，Yang［9］開發(fā)了一個三維非穩(wěn)態(tài)模型來研究微重力環(huán)境下球形玻璃容器內(nèi)液體表面的運動。目前，關(guān)于水平方管的研究結(jié)果包含了微重力與常重力環(huán)境下，氣液及液液的流動特性，本文通過數(shù)值模擬方法，建立水平方管內(nèi)的氣液兩相流三維非穩(wěn)態(tài)模型，對部分重力環(huán)境下空氣－水和制冷劑R134a 蒸汽－液體兩相的流動的壓降特性進行數(shù)值研究。

1 物理模型

根據(jù)實驗中應(yīng)用成熟的典型設(shè)計建立了數(shù)值計算的物理模型，如圖1。氣體和液體由交叉部分流入并混合，隨后流入水平主管道。設(shè)定液體沿z軸方向引入，氣體沿y 軸方向引入(微重力環(huán)境下的兩相流動不依賴于流動方向［6］)。由于當(dāng)試驗段長徑比大于50 時，可以忽略進口效應(yīng)對失重條件下氣液兩相流型的影響［10］，故計算過程重點關(guān)注于50D～60D 間管段的壓降變化。

圖1 物理模型Fig.1 Physical model

2 數(shù)值模擬

(1)控制方程

流體體積函數(shù)模型(volume of fluid，VOF)已經(jīng)被廣泛應(yīng)用于非穩(wěn)態(tài)氣液兩相流動計算［11］。本文基于VOF 模型，采用FLUENT 內(nèi)置跟蹤界面的界面重構(gòu)法(piecewise linear interface construction，PLIC)計算氣液交界面的平流流態(tài)［12］，通過連續(xù)表面張力模型(continuous surface force，CSF)計算表面張力對動量方程的影響［13］，應(yīng)用RNG k－ε湍流模型結(jié)合標準壁面函數(shù)研究湍流的影響，k 和ε 由整個兩相流域確定。研究中假定兩相流體不易混溶，不可壓縮，兩相間沒有質(zhì)量輸運，不考慮過程中的熱交換。

VOF 模型中，每個控制體內(nèi)所有相的體積分數(shù)的和為1。將氣相的體積分數(shù)記為αG，液相的體積分數(shù)記為αL，則:

質(zhì)量守恒方程和動量守恒方程:

每個控制體內(nèi)的密度和黏度分別定義為

式中:σ是表面張力;ρ 是體積平均密度;表面曲率κ 為

計算區(qū)域如圖1所示，采用均勻六面體網(wǎng)格劃分。作者對比了邊長為0.05 mm，0.1 mm，0.2 mm 和0.5 mm 四種不同尺寸網(wǎng)格的計算精度和計算成本。圖2 的對比結(jié)果顯示，與0.2 mm 和0.5 mm 網(wǎng)格相比，0.1 mm 網(wǎng)格的計算精度顯著提高;與0.1 mm 網(wǎng)格相比，0.05 mm 網(wǎng)格的計算精度提高不明顯，且需要耗費更多的計算時間。故選擇0.1 mm 網(wǎng)格，網(wǎng)格總數(shù)約為308 000 個，時間步長為10－5s。

(2)邊界條件

入口邊界條件選擇方向垂直于入口邊界的恒定速度入口，出口邊界條件是壓力出口。物面邊界條件是無滑移固壁邊界。初始時刻管道內(nèi)充滿處于靜止狀態(tài)的氣體。

(3)控制方程的離散

控制方程的對流項采用二階迎風(fēng)格式離散，壓力方程的差分格式選擇二階迎風(fēng)格式來離散控制方程的對流項，選擇交錯壓力格式(Pressure Staggering Option，PRESTO)作為壓力方程的差分格式，選擇壓力的隱式算子分割法(Pressure Implicit with Splitting of Operators，PISO)進行壓力－速度耦合計算。以殘差值全部小于10－5作為收斂標準。

圖2 網(wǎng)格尺寸分別為0.05 mm，0.1 mm，0.2 mm 和0.5 mm 時液體入口處的相分布圖Fig.2 Contours of phase at the entrance of liquid with different mesh sizes: 0.05 mm，0.1 mm，0.2 mm and 0.5 mm

3 結(jié)果與討論

(1)數(shù)值方法驗證

表1 實驗物性參數(shù)Tab.1 Physical parameters in experiments

表2 實驗參數(shù)范圍Tab.2 Parameter range in experiments

根據(jù)趙建福等人［6］實驗中的操作工況(實驗材料物性參數(shù)列于表1，實驗參數(shù)范圍列于表2)進行數(shù)值模擬。通過統(tǒng)計不同流型出現(xiàn)的工況得到了流型分布圖，數(shù)值計算結(jié)果同實驗結(jié)果的對比如圖3所示，圖中字母代表的流型列于表3。圖中的點為數(shù)值模擬得到的結(jié)果，作者在0.1g0重力環(huán)境下計算了30 組數(shù)據(jù)用于數(shù)值方法驗證。表2中，Bond 數(shù)、Weber 數(shù)和Reynold 數(shù)均為無量綱數(shù)，計算式依次為式中ρ 為密度;a 為剩余加速度;u 為表觀流速;d為特征尺度;σ 為表面張力。除與趙建福等［6］的實驗結(jié)果對比外，圖中還標志了Weisman 模型［14］和Lee 模型［15］預(yù)測的微重力環(huán)境下氣液流型的轉(zhuǎn)變界限。由對比可知，模擬結(jié)果與實驗結(jié)果和模型的預(yù)測結(jié)果都吻合較好。

圖3 數(shù)值計算結(jié)果與實驗結(jié)果對比Fig.3 Comparisons between numerical results and experimental results

表3 字母代表的流型Tab.3 Flow patterns represented by letters

(2)部分重力環(huán)境下的壓降分布

模擬選用材料的物理參數(shù)列于表4，其中空氣和水的溫度是20℃，壓力環(huán)境是0.10 MPa;R134a蒸汽和液體的溫度是飽和溫度30℃，壓力環(huán)境是0.77 MPa;重力場分別選取g0(常重力，g0=9.8 m/s2)，0.38g0(馬丁重力)和0.17g0(月球重力)進行計算。數(shù)值計算過程中氣液兩相工況參數(shù)范圍列于表5，其中，U 為表觀速度，I 為表觀動量通量(I =ρU2)，We 為Weber 數(shù)，Re 為Reynolds 數(shù)，F(xiàn)r 為Froude 數(shù)。

一般地，管道內(nèi)的流動壓降主要由沿程摩擦壓降Δpf、局部壓降Δpc、沿程加速壓降Δpa和重力壓降Δpg四部分組成［16］，即

表4 計算物性參數(shù)Tab.4 Physical parameters in numerical simulations

表5 數(shù)值計算選取的工況范圍Tab.5 Parameter range in numerical simulations

由于部分重力環(huán)境下水平管內(nèi)氣液兩相絕熱流動可忽略Δpc，Δpa和Δpg，則Δp = Δpf。

圖4所示為水平方管內(nèi)不同氣液速度下，空氣－水和R134a 蒸汽－液體兩相流在部分重力和常重力環(huán)境下的壓降分布。圖中I 為表觀動量通量，下標G 和L 依次代表氣相和液相。由圖可知，兩種混合物的壓降均隨著氣、液速度的增大而增大，這一現(xiàn)象與氣液兩相流壓降的一般規(guī)律相符。相同工況時，部分重力下的壓降大于常重力下的壓降。其原因是與常重力環(huán)境相比，部分重力環(huán)境下的浮力減小，滑移比的降低導(dǎo)致液相流速的增大，氣泡運動引起的紊動幅度減低。由于水的粘度是R134a 液體粘度的五倍，故與空氣－水兩相流相比，R134a 蒸汽－液體的壓降較小，這種現(xiàn)象隨著重力環(huán)境的減小而變得越明顯。

(3)壓降的對比及修正

本節(jié)將部分重力環(huán)境下空氣－ 水兩相流和R134a 蒸汽－液體兩相流壓降的數(shù)值結(jié)果與各模型的預(yù)測值進行對比分析。圖5所示為數(shù)值結(jié)果和均相流模型［17］的對比。由圖可知，模型的預(yù)測值整體偏大，在摩擦壓降較小時其偏差更大。這是因為，均相流模型是在常重力環(huán)境下提出的壓降預(yù)測模型，在部分重力環(huán)境下由于氣液相對運動強度減弱，使得氣液與壁面摩擦阻力顯著減小，導(dǎo)致均相流模型的預(yù)測值出現(xiàn)了較大偏差。大約61.24 %的數(shù)據(jù)高于+30 %的偏差線，均相流模型不能準確地預(yù)測部分重力環(huán)境下氣液兩相壓降的數(shù)值計算結(jié)果。

圖4 不同重力環(huán)境下空氣－水兩相流、R134a 蒸汽－液體兩相流的壓降對比Fig.4 Comparisons of pressure drop for air-water two phase flow and R134a vapor-liquid two phase flow under different gravity

圖5 模擬值與均相流模型預(yù)測值對比Fig.5 Comparison of numerical result with homogeneous flow model

數(shù)值結(jié)果和Friedel 模型［18］間的對比如圖6所示。由圖可知，F(xiàn)riedel 模型的預(yù)測值比數(shù)值結(jié)果偏大，約18.55 %的數(shù)據(jù)超出了+30 %的誤差線。Friedel 模型也難以獲得令人滿意的預(yù)測結(jié)果。

數(shù)值結(jié)果與Chisholm 模型［19］預(yù)測值的比較結(jié)果如圖7所示。相對均相流模型和Friedel 模型而言，Chisholm 模型的預(yù)測值較為理想，總體偏差分布在+30 %偏差線以內(nèi)。為了得到較理想的模型預(yù)測結(jié)果，需要對現(xiàn)有的模型進行修正。

圖6 模擬值與Friedel 模型預(yù)測值對比Fig.6 Comparison of numerical result with Friedel model

圖7 模擬值與Chisholm 模型預(yù)測值對比Fig.7 Comparison of numerical result with Chisholm model

為方便應(yīng)用，Chisholm 將Lockhart-Martinelli關(guān)聯(lián)式簡寫為

式中:關(guān)聯(lián)參數(shù)C 的取值取決于氣、液流動狀態(tài)，具體如表6所示。

表6 Tabular 常數(shù)的取值Tab.6 Values of tabular constant

圖8 為數(shù)值計算結(jié)果與該關(guān)聯(lián)式預(yù)測值的對比。由圖可知，氣、液兩相流動均為層流流態(tài)時，約36.11 %的數(shù)值結(jié)果超出了關(guān)聯(lián)式的預(yù)測值;氣、液兩相流動均為湍流流態(tài)時，約40.52 %的數(shù)值結(jié)果超出了關(guān)聯(lián)式的預(yù)測值。因此，根據(jù)氣液兩相Reynolds 數(shù)計算C 值難以對部分重力環(huán)境下結(jié)果進行準確預(yù)測。大量文獻研究中也顯示，根據(jù)氣液兩相Reynolds 數(shù)對關(guān)聯(lián)參數(shù)C 進行取值，難以全面準確地體現(xiàn)流動的實際情況和參數(shù)，故需要對現(xiàn)有模型進行修正。

圖8 數(shù)值結(jié)果與Chisholm 的分區(qū)對比Fig.8 Comparison of numerical result with Chisholm correlation

按照分液相Reynolds 數(shù)對流動重新分區(qū):層流區(qū)為ReL＜1 200;過渡區(qū)為1 200 ＜ReL＜2 100;湍流區(qū)為ReL＞2 100。由于Chisholm 關(guān)聯(lián)式中兩相相互作用項中X 指數(shù)并不為1［20］，故將該關(guān)聯(lián)式修正為

對數(shù)值結(jié)果在層流區(qū)、過渡區(qū)和湍流區(qū)的分別進行最小二乘擬合，得到m =1.31，得到的關(guān)系式列于表7。

表7 Tabular 常數(shù)的取值Tab.7 Values of tabular constant

圖9 根據(jù)液相Reynolds 數(shù)分區(qū)時ΦL2與X 的關(guān)系Fig.9 Relationship between ΦL2 and X divided according to liquid Reynolds number

圖10 模擬值與Chisholm 模型預(yù)測值對比Fig.10 Comparison of numerical result with Chisholm model

圖9 為根據(jù)分液相Reynolds 數(shù)對流動進行重新分區(qū)的結(jié)果。根據(jù)修正后的Chisholm 關(guān)聯(lián)式，有89.7 %的數(shù)值結(jié)果分布在模型預(yù)測范圍內(nèi)。數(shù)值計算結(jié)果與擬合公式預(yù)測值的比較如圖10所示，有93.1 %的點落在±25 %的誤差帶以內(nèi)，與數(shù)值計算結(jié)果有較高的吻合度。

4 結(jié)論

(1)與常重力環(huán)境相比，部分重力環(huán)境下的氣液兩相流壓降較大;與常重力相同，部分重力環(huán)境下兩相流的壓降隨氣、液入口速度的提高而增大;相同工況下，由于R134a 的粘度僅為水的1/5，故R134a 蒸汽－液體的壓降明顯低于空氣－水兩相流的壓降。

(2)對于部分重力環(huán)境下的氣液兩相流動，均相流模型和Friedel 模型的預(yù)測值與數(shù)值計算結(jié)果存在較大偏差，Chisholm 模型的預(yù)測值偏差較小，但仍然難以得到滿意的預(yù)測結(jié)果。綜合分析數(shù)值計算結(jié)果，以分液相Reynolds 數(shù)為基礎(chǔ)重新對流動分區(qū)，同時對Chisholm 關(guān)聯(lián)式修正。對比顯示，修正后的關(guān)聯(lián)式能較好地預(yù)測部分重力環(huán)境下氣液兩相流的壓降。

［1］Rezkallah K S.Recent progress in the studies of twophase flow at microgravity conditions［J］.Advances in Space Research，1995，16:123–132.

［2］Hurlbert K M，White L C，Best F R，et al.Scaling two-phase flows to Mars and Moon gravity conditions［J］.International Journal of Multiphase Flow，2004，30:351–368.

［3］Park C，Sunada E.Vapor compression hybrid twophase loop technology for Lunar surface applications［J］.Journal of AIP Conference Proceedings，2008，969，:37–43.

［4］趙建福.微重力條件下氣/液兩相流流型的研究進展［J］.力學(xué)進展，1999，29 (3):369－382.

［5］Zhao J F，Xie J C，Lin H，et al.Experimental studies on two-phase flow patterns aboard the Mir space station［J］.International Journal of Multiphase Flow，2001，27 (11):1931－1944.

［6］趙建福，解京昌，林海等.部分重力條件下氣液兩相流型研究［J］.工程熱物理學(xué)報，2004，25(1):85－87.

［7］Clarke N N，Rezkallah K S.A study of drift velocity in bubbly two-phase flow under microgravity conditions［J］.International Journal of Multiphase Flow，2001，27:1533–1554.

［8］Gabriel K S.Microgravity two-phase flow and heat transfer［M］.Springer，Berlin，2007:.

［9］Yang A.Investigation of liquid–gas interfacial shapes in reduced gravitational environments［J］.International Journal of Mechanical Science，2008，50:1304－1315.

［10］趙建福，林海，解京昌，等.失重飛機搭載氣液兩相流實驗研究［J］.空間科學(xué)學(xué)報，2000，20(4):340－347.

［11］Hirt C W，Nichols B D.Volume of fluid (VOF)method for the dynamics of free boundaries［J］.Computer Physics，1981，46:201－225.

［12］Morton K W，Baibnes M J.Numerical methods for fluid dynamics［M］.New York:Academic Press，1982:273－274.

［13］Brackbill J U，Kothe D B，Zemach C.A continuum method for modeling Surface tension［J］.Computer Physics，1992，100:335–354.

［14］Weisman J，Kang Y.Flow pattern transitions in vertical and upwards inclined lines［J］.International Journal of Multiphase Flow，1981，7 (3):271－291.

［15］Lee J.Scaling analysis of gas-liquid two-phase flow pattern in microgravity.In:31stAerospace Sciences Meeting ＆ Exhibit，Reno，NV，1993.

［16］周云龍，孫斌，陳飛.氣液兩相流型智能識別理論及方法［M］.北京:科學(xué)出版社，2007.112－117.

［17］周云龍，郭婷婷.高等流體力學(xué)［M］.北京:中國電力出版社，2008.171－173.

［18］Friedel L.Improved friction pressure drop correlations for horizontal and vertical two-phase pipe flow［C］.Ispra，Italy:European Two-phase Group Meeting，1979.

［19］Chisholm D.Two-Phase Flow in Pipelines and Heat Exchanger［M］.New York:Longman，1983.

［20］徐建鴻，駱廣生，陳桂光，等.液－液微尺度混合體系的傳質(zhì)模型［J］.化工學(xué)報，2005，56 (3):435–440.