蒸汽發生器傳熱管一、二次側耦合換熱及管外過冷沸騰數值研究

王成龍，叢騰龍，王澤勇，田文喜，秋穗正，蘇光輝，安洪振

(1.西安交通大學 能源與動力工程學院，陜西 西安 710049；2.環境保護部 核與輻射安全中心，北京 100082)

在蒸汽發生器運行過程中，二次側常伴有復雜的兩相流動，會明顯改變工質流動和傳熱特性，對蒸汽發生器的安全性與可靠性造成很大影響。近年來，蒸汽發生器二次側兩相流動沸騰現象研究備受關注。

目前，對蒸汽發生器二次側流動沸騰的研究大多基于實驗，所得實驗關聯式的局限性較大。也有研究者采用數值模擬方法對蒸汽發生器二次側流動沸騰進行了研究，大多采用RETRAN、RELAP5和THERMIT等系統程序，但這些程序基本采用一維或簡化多維模型，并不能提供蒸汽發生器二次側的局部流場信息。鑒于此，本工作應用大型商用CFD軟件ANSYS CFX 12.0，采用兩流體歐拉數學模型并結合Kurual等[1]提出的RPI壁面沸騰模型，對蒸汽發生器二次側帶梅花孔板的兩相流動進行模擬，并與一次側、管壁進行耦合傳熱計算。

1 數學物理模型

過冷沸騰的多維數值分析計算常采用兩流體歐拉數學模型進行，其對氣相和液相分別求解一套質量、動量、能量守恒方程，同時考慮氣液兩相間的傳熱、傳質和動量交換過程。

1.1 兩流體歐拉數學模型

質量守恒方程：

ρα·(rαραUα

(1)

動量守恒方程：

ραUα·(rα(ραUα?Uα))=

(2)

能量守恒方程：

ραhα·(rα(ραUαhα-λα

(3)

1.2 相間動量交換模型

(4)

(5)

式中：db為氣泡平均直徑，m；拽力系數CD取決于氣泡雷諾數Reb，本文采用Ishii-Zuber關系式[2]計算。

(6)

1.3 相間熱量傳輸模型

相間傳熱采用雙熱阻模型，用界面傳熱系數來描述，習慣上采用無量綱努賽爾數Nu表示傳熱系數。在液相側，采用Ranz-Marshall經驗關系式[4]計算努賽爾數：

(7)

在氣相側，采用零熱阻模型，其換熱系數hα為無窮大，等效于交界面溫度，等于氣相飽和溫度。

1.4 相間質量傳輸模型和RPI壁面沸騰模型

相間質量傳輸采用熱相變模型，而在壁面上的沸騰現象采用RPI壁面沸騰模型進行模擬計算。RPI壁面沸騰模型主要包括兩部分：壁面熱流密度分配模型和輔助模型。

1) 壁面熱流密度分配模型

Judd等[5]的研究結果表明：過冷沸騰時，壁面上的熱流密度qw分為3部分：液體單相對流換熱帶走的熱流密度qc，液相蒸發產生氣泡時帶走的潛熱熱流密度qe以及淬滅熱流密度qq。

qw=qc+qe+qq

(8)

qc的計算公式為：

qc=hcA1(Tw-Tf)

(9)

式中：hc為單相對流換熱系數；A1為壁面上液相所占面積；Tw為壁面溫度；Tf為液相溫度。

qe表示為：

πdbwρgnfhfg/6

(10)

qq采用Mikic等[6]提出的關系式進行計算：

(11)

式中：A2為壁面上氣相所占面積份額，A2=1-A1；τw為氣泡等待時間；cpf為液相比定壓熱容。

2) 輔助模型

氣泡成核密度n采用Lemmert-Chawla經驗關系式[7]進行計算：

n=[210(Tw-Tsat)]1.805

(12)

式中，Tsat為相應壓力下的飽和溫度

氣泡脫離直徑dbw采用Tolubinski關系式[8]進行求解：

(13)

式中：dref為參考氣泡直徑，dref=0.6 mm；dmax為最大氣泡直徑，dmax=1.4 mm；ΔTref為參考溫差，ΔTref=45 K；ΔTsub為近壁面過冷度。

氣泡脫離頻率f采用Cole關系式[9]計算：

(14)

式中：g為重力加速度；系數CD近似等于1。

氣泡等待時間τw采用Tolubinski關系式[8]進行求解：

τw=0.8/f

(15)

氣泡平均直徑db采用Anglart線性關系式[10]進行計算：

(16)

式中，d0和d1分別為參考過冷度θ0和θ1下的參考氣泡直徑。

1.5 湍流模型

由于梅花孔板附近有一部分流動為層流，在這部分區域采用k-ε模型將會導致計算發散，因此對于液相采用適應性良好的SST模型，而氣相采用彌散相零方程模型。

2 網格劃分及邊界條件

本文所選取的計算區域為帶有支撐板的蒸汽發生器傳熱管束，包括一次側、管壁及二次側。支撐板為三葉梅花孔結構，管束排列為三角形。計算區域示意圖如圖1所示。高溫的一次側流體從底端流入，將熱量傳遞給管壁，二次側流體從底端向上流經梅花孔板，最終在一未知高度處發生過冷沸騰。網格劃分如圖2所示，網格質量大于0.5。

網格進口邊界條件為速度進口邊界條件，出口邊界條件為壓力出口。一次側壁面邊界條件設置為無滑移壁面。二次側壁面邊界條件(包括加熱壁面和梅花孔板)，對液相設置為無滑移壁面邊界條件，對氣相設置為滑移壁面邊界條件。對于管壁，為便于收斂，底面設置為二次側進口溫度，頂面設置為絕熱壁面。其他邊界設置為對稱邊界條件。一次側與壁面的網格連接方式為GGI模式，二次側與壁面的網格連接方式為1∶1模式。

a——支撐板梅花孔結構示意圖；b——傳熱管束三維結構示意圖

圖2 蒸汽發生器傳熱管束網格劃分示意圖

3 結果分析及模型驗證

3.1 模型驗證及網格敏感性分析

本文采用Bartolomej等[11]做的圓管內過冷沸騰實驗來驗證模型的適用性和準確性。實驗圓管長2 m，內徑15.4 mm，壁面熱流密度5.7×105W/m2，進口質量流量900.0 kg/(s·m2)，壓力4.5 MPa，進口過冷度58.2 K，具體幾何結構如圖3所示。計算結果與實驗值的比較如圖4所示。圖4結果表明，本模型對過冷沸騰的模擬與實驗符合較好，其最大相對誤差為3.2%，可認為RPI壁面沸騰模型能較好地模擬過冷沸騰現象。

圖3 標準題幾何模型

圖4 計算結果與實驗值比較

為驗證網格的獨立解，分別對3套網格進行了計算，計算結果列于表1。可看出，隨著網格數量的增加，二次側出口空泡份額和出口溫度趨于定值，從而驗證了網格的獨立解。本次計算中采用第2套網格(1 118 040)進行計算。

表1 網格敏感性分析

3.2 結果分析

計算工況列于表2。管壁及一、二次側工質的物性由查表所得。

表2 計算工況

圖5為二次側流體流線分布。由于二次側流體通過梅花孔板時，流通截面突縮，導致流體在此突然加速，最大流速為1.57 m/s。梅花孔處由于流通面積逐漸縮小及其獨特的結構，導致流體在流經梅花孔前后產生逆時針旋轉流動，如圖5b所示。同時還可看出，在梅花孔板窄縫區，由于流體受到兩側壁面(管壁和孔板壁面)摩擦力的作用，其流動速度非常低，為0.004 m/s。

傳熱管束對稱面處一、二次側流體溫度及其梅花孔板局部溫度分布如圖6所示。由圖6a可見，一次側高溫流體由進口流入，不斷將熱量傳入管壁，管壁將熱量傳遞給二次側流體，最終導致一次側流體溫度降低，其出口平均溫度為563.47 K。二次側流體隨熱量的不斷導入，其溫度逐漸上升，出口平均溫度為543.11 K。由圖6b可見，由于二次側流體流過梅花孔板處產生逆時針渦旋，其出口處的溫度分布不再對稱，區域A處的溫度略高，而區域B處的略低。由于一、二次側的耦合效應，一次側溫度也呈現不均勻分布，A側的換熱系數高導致一次側溫度下降較大，而B處的換熱系數低導致一次側溫度下降較少。由圖6c、d可見，梅花孔板窄縫區的流速很低，導致流體與壁面的換熱由原來的對流換熱方式轉變為以導熱方式為主的換熱，此處熱量聚集導致溫度出現局部最大，為552.37 K，超過二次側流體飽和溫度546.54 K，因此，此處很容易發生液相干涸，導致臨界熱流密度的發生。其他區域的換熱以對流換熱方式為主，靠近壁面的溫度較高，而中心區域的溫度較低，如圖6d所示。

a——二次側整體流速分布；b——二次側流速俯視圖

a——一、二次側溫度分布；b——一、二次側出口溫度分布；c——梅花孔板側面溫度分布；d——梅花孔板y=0截面處溫度分布

蒸汽發生器二次側對稱面及梅花孔板局部空泡份額分布示于圖7。由圖7a可見，二次側流體在離進口約0.23 m處發生過冷沸騰，隨著熱量的不斷導入，截面平均空泡份額逐漸增大，而壁面溫度卻增加緩慢，約為552.12 K。由于流體逆時針旋轉，導致壁面熱流密度分布不均，從而造成截面空泡份額分布不均(圖7b)。由圖7b可見，A區域的平均空泡份額明顯大于B區域的，出口平均空泡份額為0.091。由圖7c可見，流體流經梅花孔板前時近壁面處空泡份額明顯大于流經梅花孔板之后，這是由于梅花孔的結構使液相和氣相發生強烈攪混，使得氣泡在過冷液相中冷凝的概率變大，最終導致近壁面空泡份額下降。而在孔板窄縫區，由于此處液體的流速極慢，導致此處液相干涸，空泡份額存在一局部峰值，約為0.19。圖7d～f示出流體流過梅花孔板不同位置截面處的空泡份額分布，可看出，流體從流進梅花孔板到流出梅花孔板，其截面平均空泡份額依次增大。在截面y=7 mm處存在梅花孔板窄縫區域，存在空泡份額峰值。

a——二次側空泡份額分布；b——二次側出口處空泡份額分布；c——梅花孔板處空泡份額分布；d——梅花孔板y=-7 mm處空泡份額分布；e——梅花孔板y=0 mm處空泡份額分布；f——梅花孔板y=7 mm處空泡份額分布

圖8 不同高度處截面平均空泡份額分布

圖8示出不同高度處截面平均空泡份額分布。由圖8可看出，流體流經梅花孔板之前，截面平均空泡份額很低。這時主要為單相液體流動，其壁面上的換熱以對流換熱為主，而淬滅換熱和沸騰換熱所占份額很小。流體剛流入梅花孔板，空泡份額降低，這是因為此處液相和氣相發生了強烈攪混。隨后，空泡份額迅速上升達到0.023，而在孔板出口處空泡份額再次突降，這由流通面積突擴及出口處氣液兩相的強力攪混兩方面因素造成。之后，由于熱量不斷導入，截面平均空泡份額迅速增加，達到0.091。

圖9 流體壁面平均換熱系數隨高度的變化

圖9示出一、二次側流體壁面平均換熱系數隨高度的變化。由圖9可看出，一次側的換熱系數隨著高度的增加而增大，在出口處達到最大(50 843.90 W/(m2·K))，梅花孔板的存在并未對一次側換熱系數造成影響。二次側換熱系數的變化較為復雜，流體換熱系數剛開始隨著高度的增加而減小，這主要是由于此時流體流動主要為單相液體流動，進口處溫差(主流與壁面)較大，從而導致較高的壁面熱流密度，造成較大的換熱系數。流體流過梅花孔板處，傳熱系數突然上升，達到最大值(27 057.70 W/(m2·K))。這主要是由兩方面因素造成的：一是截面積的突縮導致流速突然增加，使得對流換熱效應大幅增強；二是由于此處的過冷沸騰相對劇烈，增強了換熱。此后，壁面換熱系數隨著高度的增加先有所減小，后逐漸增大。

表3列出蒸汽發生器一、二次側耦合傳熱過冷沸騰關鍵參數的分析結果。

表3 過冷沸騰的分析結果

4 結論

1) 采用RPI壁面沸騰模型對蒸汽發生器二次側進行過冷沸騰模擬，準確預測了沸騰起始點的位置，在離進口約0.23 m處。利用該模型計算的結果與Bartolomej等[11]的實驗結果進行對比，符合良好。

2) 梅花孔板的存在使二次側流體發生了逆時針渦旋，最終導致一次側、二次側的流場、溫度場和空泡份額分布不均勻。

3) 梅花孔板窄縫區使得該處的流速極小，導致液相干涸，空泡份額在此處有一局部峰值，易導致傳熱惡化，管束損毀。因此，利用CFD對梅花孔板區域進行數值模擬對梅花孔板結構設計具有重要意義。

參考文獻：

[1] KURUAL N, PODOWSKI M Z. On the modeling of multidimensional effects in boiling channels[C]∥ANS Proceedings of 27th National Heat Transfer Conference. Minneapolis: [s. n.], 1991.

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

[3] BURNS A D, FRANK T, HAMILL I, et al. The Favre averaged drag model for turbulent dispersion in Eulerian multi-phase flows[C]∥5th International Conference on Multiphase Flow, ICMF-2004. Yokohama, Japan: [s. n.], 2004.

[4] RANZ W E, MARSHALL W R. Evaporation form drops[J]. Chem Eng Prog, 1952, 48(3): 141-146.

[5] JUDD R L, HWANG K S. A comprehensive model for nucleate pool boiling heat transfer including microlayer evaporation[J]. ASME J Heat Transfer, 1976, 94(4): 623-629.

[6] MIKIC B, ROHSENOW W M. A new correlation of pool-boiling data including the fact of heating surface characteristics[J]. ASME J Heat Transfer, 1969, 91(2): 245-250.

[7] LEMMERT M, CHAWLA J M. Influence of flow velocity on surface boiling heat transfer coefficient[M]∥Heat Transfer and Boiling. American: Academic Press, 1977.

[8] TOLUBINSKI V I, KOSTANCHUK D M. Vapor bubbles growth rate and heat transfer intensity at subcooled water boiling[C]∥4th International Heat Transfer Conference. Paris, France: [s. n.], 1970.

[9] COLE R. A photographic study of pool boiling in the region of CHF[J]. AIChE J, 1960, 6(4): 533-542.

[10] ANGLART H, NYLUND O. CFD application to prediction of void distribution in two-phase bubbly flows in rod boundles[J]. Nucl Eng Des, 1996, 163(1-2): 81-98.

[11] BARTOLOMEJ G, BRANTOV V, MOLOCHNIKOV Y S, et al. An experimental investigation of true volumetric vapour content with subcooled boiling in tubes[J]. Thermal Engineering, 1982, 29(3): 132-135.