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

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

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

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

1 數(shù)學物理模型

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

1.1 兩流體歐拉數(shù)學模型

質(zhì)量守恒方程：

ρα·(rαραUα

(1)

動量守恒方程：

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

(2)

能量守恒方程：

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

(3)

1.2 相間動量交換模型

(4)

(5)

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

(6)

1.3 相間熱量傳輸模型

相間傳熱采用雙熱阻模型，用界面?zhèn)鳠嵯禂?shù)來描述，習慣上采用無量綱努賽爾數(shù)Nu表示傳熱系數(shù)。在液相側(cè)，采用Ranz-Marshall經(jīng)驗關(guān)系式[4]計算努賽爾數(shù)：

(7)

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

1.4 相間質(zhì)量傳輸模型和RPI壁面沸騰模型

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

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

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

qw=qc+qe+qq

(8)

qc的計算公式為：

qc=hcA1(Tw-Tf)

(9)

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

qe表示為：

πdbwρgnfhfg/6

(10)

qq采用Mikic等[6]提出的關(guān)系式進行計算：

(11)

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

2) 輔助模型

氣泡成核密度n采用Lemmert-Chawla經(jīng)驗關(guān)系式[7]進行計算：

n=[210(Tw-Tsat)]1.805

(12)

式中，Tsat為相應(yīng)壓力下的飽和溫度

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

(13)

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

氣泡脫離頻率f采用Cole關(guān)系式[9]計算：

(14)

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

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

τw=0.8/f

(15)

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

(16)

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

1.5 湍流模型

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

2 網(wǎng)格劃分及邊界條件

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

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

a——支撐板梅花孔結(jié)構(gòu)示意圖；b——傳熱管束三維結(jié)構(gòu)示意圖

圖2 蒸汽發(fā)生器傳熱管束網(wǎng)格劃分示意圖

3 結(jié)果分析及模型驗證

3.1 模型驗證及網(wǎng)格敏感性分析

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

圖3 標準題幾何模型

圖4 計算結(jié)果與實驗值比較

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

表1 網(wǎng)格敏感性分析

3.2 結(jié)果分析

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

表2 計算工況

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

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

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

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

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

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

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

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

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

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

表3列出蒸汽發(fā)生器一、二次側(cè)耦合傳熱過冷沸騰關(guān)鍵參數(shù)的分析結(jié)果。

表3 過冷沸騰的分析結(jié)果

4 結(jié)論

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

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

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

參考文獻：

[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.