基于非平衡熱力學模型的管內沸騰過程模擬

馮留海， 聶永廣， 王江云， 雷曉東， 王　娟， 毛　羽

(1.中國石油大學 重質油國家重點實驗室， 北京 102249； 2.北京低碳清潔能源研究所， 北京 102209；3.新奧科技發展有限公司， 河北 廊坊 065001； 4.中機康元糧油裝備(北京)有限公司， 北京 100083)

基于非平衡熱力學模型的管內沸騰過程模擬

馮留海1,2， 聶永廣1,3， 王江云1， 雷曉東4， 王娟1， 毛羽1

(1.中國石油大學 重質油國家重點實驗室， 北京 102249； 2.北京低碳清潔能源研究所， 北京 102209；3.新奧科技發展有限公司， 河北 廊坊 065001； 4.中機康元糧油裝備(北京)有限公司， 北京 100083)

摘要：為了更準確地模擬壁面沸騰過程，根據Lavievile非平衡熱力學模型對沸騰流動的壁面熱流率進行了劃分，通過調節函數結合Sato模型修正了沸騰流動中的混合k-ε模型，建立了同時適用于管內過冷沸騰與整體沸騰流動階段的曳力模型。使用上述模型對環狀豎直管內沸騰流動過程進行數值模擬，計算結果與實驗值吻合，可以用來模擬整個管內沸騰流動過程。數值模擬顯示，沸騰流動可以明顯增大管壁的傳熱系數，但當管壁處蒸氣體積分數超過0.6時，會出現傳熱惡化現象。

關鍵詞：沸騰流動； 非平衡熱力學； 曳力模型； 傳熱系數； 數值模擬

管內沸騰流動是石油、化工、能源動力等工業領域的一種常見現象，例如石油化工行業常用的管式加熱爐和管殼式換熱器中的管內沸騰流動。沸騰流動可以顯著提高管內傳熱效率，所以受到人們的廣泛關注。早在20世紀60年代，Bartolomei等[1]就實驗研究了空泡分率對垂直管內高壓水沸騰過程的影響。隨著測量技術的發展，Bartel[2]、Lee等[3]對各種工況下的管內沸騰過程進行了更加詳細的測量。一方面希望實驗獲得管內沸騰過程中詳細的兩相流局部參數特性，另一方面希望通過實驗測量結合理論方法歸納出氣泡脫離直徑、成核頻率和氣泡停留時間等管內沸騰流動中的關鍵參數，從而完善數值模型，使其可以準確模擬管內沸騰過程。Kurul等[4]最早基于雙流體模型對壁面熱流率進行了劃分，成功建立了RPI模型。因其考慮了壁面熱流率的分配情況，所以能更真實地反映管內高壓水過冷沸騰過程，并在相關領域得到廣泛采用。此后，大量學者不斷完善RPI的子模型。Anglart等[5]采用RPI模型研究了核燃料分配器中的過冷沸騰流動過程，Kocar等[6-7]建立了低壓工況下水的過冷沸騰模型，李祥東等[8-9]詳細總結了面向液氮過冷沸騰流動的數值模型，竇從從等[10-11]通過數值模擬確定了垂直管內高壓高過冷水過冷沸騰流動的子模型。研究表明，RPI模型能夠較好地預測過冷沸騰流動特性，但是由于其假設蒸氣均為飽和溫度，從而影響了管內沸騰流動過程中整體沸騰階段的模擬精度。筆者在上述基礎上，并根據Lavievile的非平衡熱力學方法重新劃分了壁面熱流率[12]，通過數值計算修正了相關湍流模型和相間動量傳輸模型，最后將相關模型通過UDF植入通用CFD計算軟件Fluent，以期更好地模擬管內沸騰流動過程。

1管內沸騰過程的數學模型

1.1雙流體模型

雙流體模型假設氣、液兩相均為連續介質，根據牛頓力學法則建立各自的控制方程，兩相之間的相互作用通過控制方程中的相間傳輸相描述。在管內沸騰模型中，將液相視為連續相，蒸氣相視為擬連續相。雙流體模型是整個管內沸騰數值模型的骨架，管內沸騰流動數值模擬中所有子模型旨在封閉雙流體模型的控制方程組。適用于管內沸騰流動模擬的雙流體模型的基本控制方程組為式(1)～(3)。

質量守恒方程：

(1)

動量守恒方程：

(2)

能量守恒方程：

(3)

式中所有與相變相關的項都需要采用沸騰模型進行計算。

1.2壁面沸騰模型

管內沸騰是由加熱壁面加熱液體產生的沸騰現象，包括RPI模型在內的適用于這一過程的模型總稱為壁面沸騰模型。在這一模型中，管內的傳質傳熱主要分成2部分，即壁面處的傳質傳熱和主流區域的傳質傳熱，其中計算壁面區域傳質和傳熱的關鍵在于合理劃分壁面熱流率。

(1)壁面熱流率劃分

前人數值模擬研究采用的RPI模型大多將壁面熱流率劃分為3部分，即液相對流熱流率Ql、激冷熱流率Qq和蒸發熱流率Qe。需要說明的是，Qq是由于加熱產生的氣泡脫離壁面時補充的冷液體被瞬間加熱升溫所產生的熱流率。具體的表達式為式(4)～(6)。

Ql=hl(Tw-Tl)(1-Aq)

(4)

(5)

(6)

現有的RPI模型假設流場中的蒸氣均處于飽和溫度狀態，不考慮蒸氣在壁面處的對流換熱。本研究采用的Lavievile壁面熱流率劃分方法是在RPI模型的基礎上加入了蒸氣對流熱流率Qv，其表達式為式(7)。壁面熱流率的劃分和分配函數f(αl) 如式(8)和(9)所示。式(6)中的氣泡脫離直徑dbw使用Unal模型計算[13]，而主流區的dbw根據過冷度[8,10-11]確定。

Qv=hv(Tw-Tv)

(7)

Qw=f(αl)(Ql+Qq+Qe)+(1-f(αl))Qv

(8)

(9)

(2)相間傳質和傳熱

相間質量傳遞分為2部分。一是壁面附近液體在加熱作用下沸騰相變為蒸氣，二是蒸氣在過冷度較高的主流區域冷凝為液體。蒸氣相的傳質項表達式為式(10)，液相與蒸氣相的相間傳熱項Qlv表達式為式(11)。

(10)

Qlv=hlv(Tl-Tv)Alv

(11)

式(10)中，hlv根據Ranz-Marshall模型計算。

與RPI模型不同，因為考慮蒸氣處于非飽和溫度狀態，因此既存在液相向氣相的熱量傳遞，也存在氣相向液相的熱量傳遞。在壁面附近，壁面熱流率中的對流熱流率Ql和激冷熱流率Qq最終加入到液相能量源項，蒸發熱流率Qe和蒸氣對流熱流率Qv加入到氣相能量源項。

(3)湍流模型修正

RPI模型模擬過冷沸騰過程時，通常假設液相為湍流流動，氣相為層流流動。但隨著沸騰流動過程的發展，當管內蒸氣體積分數超過0.3時，這一假設明顯欠合理。因此，筆者使用混合k-ε模型計算兩相的混合湍流動能km和混合湍流耗散率εm，然后采用調節函數fα計算最終蒸氣湍流黏度。調節函數fα的表達式為式(12)。

(12)

式(12)中，a決定了函數曲線的變化幅度，b決定了函數曲線的閾值。圖1為調整函數fα曲線。從圖1可以看出，取a=20、b=0.3比較合適。最終有效黏度表達式為式(13)。

(13)

考慮到氣泡擾動對液相湍流的影響，使用Sato模型[14]修正液相有效湍流黏度，表達式為式(14)～(15)。

(14)

(15)

(4)相間曳力模型修正

在對管內過冷沸騰流動數值模擬中，Ishii模型[15]表現出了良好的適用性。Ishii模型原理上僅適用于氣相體積分數小于0.5的沸騰流動過程，而實際過程中隨著蒸氣量的增加，管內氣液兩相流的流型則從泡狀流逐漸過渡到霧狀流。所以，有必要建立同時適用于泡狀流動和霧狀流動的曳力模型，以研究整個管內沸騰流動過程。采用式(12)中的調節函數，將系數a、b的值分別取為20和0.5，修正的曳力系數的表達式如式(16)所示。

圖1　調整函數fα曲線

K=fα·Kbubble+(1-fα)·Kdroplet

(16)

式(16)中的Kbubble采用Ishii模型計算，而Kdroplet則采用適用性較好的Schiller模型計算。

2加熱管幾何模型及模擬工況

加熱管數值計算的幾何模型如圖2所示。模擬對象為Bartel的管內沸騰過程[2]，其內層為加熱棒，外層為玻璃管。加熱棒直徑為9.55 mm，加熱部分的長度為1500 mm(淺灰色部分)，玻璃管內徑為19.05 mm。操作工況為常壓、入口溫度367.9 K，工質為水，質量流率為700 kg/(m2·s)，壁面熱流率為188 kW/m2，出口設置為壓力出口。因為計算域為軸對稱結構，計算在柱坐標系中進行，網格劃分采用四邊形網格。壓力與速度耦合采用PCSIMPLE(Phase couple SIMPLE)算法，對流項采用QUICK差分格式。

圖2　加熱管幾何結構及尺寸

3結果與討論

圖3為Sato模型中不同Cμb取值時加熱管z=1500 mm處液相湍流黏度沿徑向的分布。由圖3看到，Cμb決定了在湍流效應中是液體剪切造成的湍流還是氣泡誘導產生的湍流起主導地位。當Cμb取值較小時，液體剪切造成的湍流占主導地位，這時液體湍流黏度在加熱壁面與管壁的中間位置最高，而在管壁附近則較低；壁面沸騰時管壁附近加熱產生的氣泡會不斷脫離壁面，破壞邊界層，增加湍流擾動，促使壁面附近的氣、液兩相快速混合，因此液相湍流黏度應該在加熱壁面附近較高。而Cμb最適宜的取值需要進一步分析其對管內溫度場和蒸氣體積分數分布造成的影響而定。

圖3　不同Cμb取值時加熱管內z=1500 mm處液

圖4為不同Cμb時 加熱管z=1500 mm處的液體溫度的徑向分布。從圖4可以看出，Cμb取值越大，表明壁面附近液體的溫度越低，說明氣泡誘導的湍流作用增強了壁面附近的熱傳遞，從而有更多的熱量被帶到主流區域。在主流區域，由于氣相體積分數的減小，氣泡產生的擾動作用減小，所以液體剪切作用產生的湍流重新占據主導地位。液體溫度會直接影響加熱壁面處氣泡脫離直徑的大小，進而影響到蒸發熱流率的大小，因此液相湍流修正對準確預測加熱管內蒸氣體積分數起到了非常重要的作用。

圖5為不同Cμb時加熱管z=1500 mm處的蒸氣體積分數沿徑向的分布。從圖5可以看出，Cμb為0.6時，計算得到的徑向蒸氣體積分數分布與實驗值最為吻合，說明Cμb=0.6時計算得到的氣泡誘導湍流作用最能反映實際情況。此外，加熱管內蒸氣體積分數最高處并非處于管壁上，而是在臨近加熱壁面的位置。這是因為在壁面潤滑力的作用下，加熱壁面附近生成的氣泡會被推離壁面，最終在壁面附近形成蒸氣富集區。

圖4　不同Cμb取值時加熱管內z=1500 mm

圖5　不同Cμb時加熱管內z=1500 mm處蒸氣體積分數

圖6為使用修正后的曳力模型和Ishii模型計算得到的加熱管內蒸氣體積分數軸向分布的對比。從圖6可以看出，在過冷沸騰段，使用修正后的曳力模型計算得到的體積分數軸向分布與Ishii模型的計算結果基本吻合，而在蒸氣體積分數超過0.4后(過冷沸騰向整體沸騰過渡)，修正后的曳力模型計算得到的體積分數高于Ishii模型的計算結果。修正后的曳力模型在管內沸騰的高體積分數區域與實驗值吻合更好。

圖6　不同曳力模型模擬的加熱管內蒸氣體積分數

圖7是加熱管壁面傳熱系數與壁面附近蒸氣體積分數軸向分布。從圖7可以看出，在加熱段開始階段，壁面傳熱系數隨著蒸氣體積分數的增加迅速增加；在蒸氣體積分數略大于0.1時，壁面傳熱系數達到峰值，之后存在一個迅速下降階段和一個穩定發展階段，隨后開始呈現線性遞減趨勢。分析認為，在管內沸騰的初始階段，壁面產生的氣泡破壞了邊界層，增強了壁面附近的流體湍流作用，導致壁面熱傳遞系數迅速升高；當壁面沸騰進一步發展時，氣泡產生的湍流也增強了氣泡和液體之間的熱傳遞和質量傳遞，而這時壁面附近液體的過冷度依然較高，使得壁面上氣泡的產生和壁面附近氣泡的冷凝達到短暫的平衡狀態，壁面附近蒸氣體積分數增加趨勢放緩，同時這種平衡狀態也弱化了壁面附近的對流熱傳遞，使得壁面傳熱系數減小；管內沸騰流動繼續發展，管壁附近液體被冷凝的氣泡充分加熱使得液體過冷度明顯降低，這時壁面附近的氣泡將不會被液體完全冷凝，加熱管壁附近的蒸氣體積分數開始持續增加，而加熱管壁附近不斷增加的蒸氣會使得壁面傳熱系數降低。從圖7可以推斷，當加熱管壁面附近的蒸氣體積分數超過0.6時，壁面傳熱系數甚至會低于單相液體加熱時的壁面傳熱系數。

圖7　加熱管壁面傳熱系數和壁面附近

4結論

(1)根據Lavievile的非平衡熱力學方法對加熱管壁面沸騰中的壁面熱流率進行劃分，使其可以適用于整個管內沸騰過程的數值計算。同時，對壁面沸騰中的兩相湍流模型修正，建立了同時適用于過冷沸騰和整體沸騰的雙流體相間曳力模型。數值模擬顯示，在使用混合k-ε模型時，Sato模型的修正系數Cμb取0.6時，計算結果與實驗值吻合較好。

(2)在豎直管內沸騰中，壁面潤滑力的作用使得蒸氣體積分數最大值出現在臨近加熱壁面的位置。沸騰現象對管壁傳熱系數有明顯的影響。在壁面處的蒸氣體積分數介于0.1～0.3之間時，壁面的傳熱系數達到最大；如果壁面處蒸氣體積分數繼續增大，壁面傳熱系數會隨之減小；當壁面處蒸氣體積分數超過0.6時，壁面傳熱系數有可能會低于單相對流傳熱的傳熱系數。

符號說明：

Aq——氣泡占據的壁面面積，m2；

Alv——相間傳熱面積，m2；

a,b——系數；

cp——液相比容，W/(m·K)；

Cμb——常系數；

db——氣泡直徑，m；

dbw——氣泡脫離直徑，m；

f——氣泡脫離頻率，Hz；

fα——調節函數；

Flv——兩相間相互作用力，包括曳力、升力、湍流耗散力和壁面潤滑力，N；

h——焓，J；

hfv——蒸發潛熱，J/kg；

hl——液相對流傳熱系數，W/(m2·K)；

hlv——相間傳熱系數，W/(m2·K)；

hv——液相傳熱系數，W/(m2·K)；

hw——壁面傳熱系數，W/(m2·K)；

k——導熱系數，W/(m·K)；

km——混合湍流強度，m2/s2；

K——曳力系數；

Kbubble——氣泡流曳力交換系數；

Kdroplet——霧狀流曳力交換系數；

mlv——液相向氣相傳質，代表液相的沸騰和蒸氣相的冷凝，kg；

Na——氣泡成核密度，m-2；

p——壓力，Pa；

Qe——蒸發熱流率，W/m2；

Ql——對流傳熱熱流率，W/m2；

Qlv——相間傳熱速率，是非相變造成的相間能量傳遞，W/m3；

Qq——激冷傳熱熱流率，W/m2；

Qw——壁面總熱流率，W/m2；

Qv——蒸氣單相對流傳熱熱流率，W/m2；

r——徑向位置，m；

R0——加熱管半徑，m；

R1——玻璃管內徑，m；

t——時間，s；

T——溫度，K；

Tsat——飽和溫度，K;

u——速度，m/s；

z——軸向，m；

α——體積分數，%；

εm——混合湍流耗散率，m2/s2；

μ——黏性系數，Pa·s；

μeff——有效湍流黏性系數，Pa·s；

μv——氣相黏性系數，Pa·s；

ρ——流體密度，kg/m3；

下標:

f，w——界面，壁面；

l，v——液相，蒸氣相。

參考文獻

[1] BARTOLOMEI G, CHANTURIYA V. Experimental study of true void fraction when boiling subcooled water in vertical tubes[J].Thermal Engineering, 1967, 14(2):123-128.

[2] BARTEL M D. Experimental investigation of subcooled boiling[D].West Lafayette: Purdue University, 1999.

[3] LEE T H, PARK G C, LEE D J. Local flow characteristics of subcooled boiling flow of water in a vertical concentric annulus[J].International Journal of Multiphase Flow, 2002, 28(8):1351-1368.

[4] KURUL N, PODOWSKI M Z. Multidimensional effects in forced convection subcooled boiling[C]//Jerusalem： Proceedings of the Ninth International Heat Transfer Conference, 1990:19-24.

[5] ANGLART H, NYLUND O, KURUL N, et al. CFD prediction of flow and phase distribution in fuel assemblies with spacers[J].Nuclear Engineering and Design, 1997, 177(1):215-228.

[6] KONCAR B, MAVKO B. CFD simulation of subcooled flow boiling at low pressure[C]//Ljubljana：Proceedings of International Conference Nuclear Energy in Central Europe, 2001:10-13.

[7] KONCAR B, KLJENAK I, MAVKO B. Modelling of local two-phase flow parameters in upward subcooled flow boiling at low pressure[J].International Journal of Heat and Mass Transfer, 2004, 47(6):1499-1513.

[8] 李祥東， 汪榮順， 黃榮國， 等. 垂直圓管內液氮流動沸騰的理論模型及數值模擬[J].化工學報， 2006, 57(3):491-497.(LI Xiangdong, WANG Rongshun, HUANG Rongguo, et al. Modelling and numerical simulation of boiling flow of liquid nitrogen in vertical tube[J].Journal of Chemical Industry and Engineering (China), 2006, 57(3):491-497.)

[9] LI Xiangdong, WEI Wei, WANG Rongshun, et al. Numerical and experimental investigation of heat transfer on heating surface during subcooled boiling flow of liquid nitrogen[J].International Journal of Heat and Mass Transfer, 2009, 52(5):1510-1516.

[10] 竇從從， 毛羽， 王娟， 等. 高壓高過冷度下過冷流動沸騰數值模擬[J].化工學報， 2010, 61(3):580-586.(DOU Congcong, MAO Yu, WANG Juan, et al. Numerical simulation of subcooled boiling flow under high pressure and high subcooling condition[J].Journal of Chemical Industry and Engineering (China), 2010, 61(3):580-586.)

[11] 竇從從， 毛羽， 王娟， 等. 高過冷度下直管內氣液兩相傳熱特性[J].化工學報， 2010, 61(9):2314-2319.(DOU Congcong, MAO Yu, WANG Juan, et al. Gas-liquid two phase flow heat transfer at high subcooling in vertical tube[J].Journal of Chemical Industry and Engineering (China), 2010, 61(9):2314-2319.)

[12] MIMOUNI S, BOUCKER M, LAVIEVILLE J, et al. Modelling and computation of cavitation and boiling bubbly flows with the NEPTUNE_CFD code[J].Nuclear Engineering and Design, 2008, 238(3):680-692.

[13] UNAL H C. Maximum bubble diameter, maximum bubble growth time and bubble growth rate[J].International Journal of Heat Mass Transfer, 1976, 19(6):643-649.

[14] SATO Y, SADATOMI M, SEKOGUCHI K. Momentum and heat transfer in two-phase bubble flow—I Theory[J].International Journal of Multiphase Flow, 1981, 7(2):167-177.

[15] ISHII M, ZUBER N. Drag coefficient and relative velocity in bubbly, droplet or particulate flows[J].AIChE Journal, 1979, 25(5):843-855.

Simulation of Boiling Flow in Vertical Tube Based on Non-Equilibrium Thermal Model

FENG Liuhai1,2, NIE Yongguang1,3, WANG Jiangyun1, LEI Xiaodong4, WANG Juan1, MAO Yu1

(1.StateKeyLaboratoryofHeavyOilProcessing,ChinaUniversityofPetroleum,Beijing102249,China;2.NationalInstituteofClean-and-LowCarbonEnergy,Beijing102209,China;3.XinAoScienceandTechnologyDevelopmentCo.Ltd.,Langfang065001,China;4.ChinaMachineryKangyuanCerealsandOilEquipment(Beijing)Co.,Ltd.,Beijing100083,China)

Abstract：To simulate the boiling process, the wall heat flux was divided based on Lavievile’s thermal non-equilibrium heat flux partition algorithm. The mixture k-ε turbulent model was improved by combination of Sato model and a modifier function, by which a new drag model was presented to simulate the full boiling flow in tube. The simulated result of boiling flow in vertical circular tube agreed with experiment data well, indicating that the model was valid in simulation of boiling flow in whole tube. It is found from the simulation that the wall heat transfer coefficient was evidently enhanced in boiling, however, heat transfer deterioration appeared when the volume fraction of vapor near the wall reached about 0.6.

Key words：boiling flow; non-equilibrium thermal dynamics; drag force model; heat transfer coefficient; numerical simulation

收稿日期：2015-04-20

基金項目：國家重點基礎研究發展計劃“973”項目(2010CB226902)和中國石油大學(北京)科研基金項目(2462015YQ0303)資助

文章編號：1001-8719(2016)03-0591-06

中圖分類號：TK124

文獻標識碼：A

doi:10.3969/j.issn.1001-8719.2016.03.021

第一作者： 馮留海，男，博士研究生，從事多相流動的數值模擬與實驗方面的研究

通訊聯系人： 王娟，女，副教授，博士，從事多相流動與燃燒過程的數值模擬與實驗研究；Tel：010-89733293；E-mail：wangjuan@cup.edu.cn；毛羽，男，教授，博士，從事多相流動及燃燒、氣固分離及液體霧化技術、化工過程裝備優化等方面的研究；Tel：010-89733293；E-mail：maoyu@cup.edu.cn