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

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

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

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

1 數(shù)學(xué)物理模型

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

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

質(zhì)量守恒方程：

ρα·(rαραUα

(1)

動(dòng)量守恒方程：

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

(2)

能量守恒方程：

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

(3)

1.2 相間動(dòng)量交換模型

(4)

(5)

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

(6)

1.3 相間熱量傳輸模型

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

(7)

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

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

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

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

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

qw=qc+qe+qq

(8)

qc的計(jì)算公式為：

qc=hcA1(Tw-Tf)

(9)

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

qe表示為：

πdbwρgnfhfg/6

(10)

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

(11)

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

2) 輔助模型

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

n=[210(Tw-Tsat)]1.805

(12)

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

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

(13)

式中：dref為參考?xì)馀葜睆?，dref=0.6 mm；dmax為最大氣泡直徑，dmax=1.4 mm；ΔTref為參考溫差，ΔTref=45 K；ΔTsub為近壁面過(guò)冷度。

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

(14)

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

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

τw=0.8/f

(15)

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

(16)

式中，d0和d1分別為參考過(guò)冷度θ0和θ1下的參考?xì)馀葜睆健?/p>

1.5 湍流模型

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

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

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

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

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

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

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

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

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

圖3 標(biāo)準(zhǔn)題幾何模型

圖4 計(jì)算結(jié)果與實(shí)驗(yàn)值比較

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

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

3.2 結(jié)果分析

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

表2 計(jì)算工況

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

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

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

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

蒸汽發(fā)生器二次側(cè)對(duì)稱面及梅花孔板局部空泡份額分布示于圖7。由圖7a可見，二次側(cè)流體在離進(jìn)口約0.23 m處發(fā)生過(guò)冷沸騰，隨著熱量的不斷導(dǎo)入，截面平均空泡份額逐漸增大，而壁面溫度卻增加緩慢，約為552.12 K。由于流體逆時(shí)針旋轉(zhuǎn)，導(dǎo)致壁面熱流密度分布不均，從而造成截面空泡份額分布不均(圖7b)。由圖7b可見，A區(qū)域的平均空泡份額明顯大于B區(qū)域的，出口平均空泡份額為0.091。由圖7c可見，流體流經(jīng)梅花孔板前時(shí)近壁面處空泡份額明顯大于流經(jīng)梅花孔板之后，這是由于梅花孔的結(jié)構(gòu)使液相和氣相發(fā)生強(qiáng)烈攪混，使得氣泡在過(guò)冷液相中冷凝的概率變大，最終導(dǎo)致近壁面空泡份額下降。而在孔板窄縫區(qū)，由于此處液體的流速極慢，導(dǎo)致此處液相干涸，空泡份額存在一局部峰值，約為0.19。圖7d～f示出流體流過(guò)梅花孔板不同位置截面處的空泡份額分布，可看出，流體從流進(jìn)梅花孔板到流出梅花孔板，其截面平均空泡份額依次增大。在截面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)梅花孔板之前，截面平均空泡份額很低。這時(shí)主要為單相液體流動(dòng)，其壁面上的換熱以對(duì)流換熱為主，而淬滅換熱和沸騰換熱所占份額很小。流體剛流入梅花孔板，空泡份額降低，這是因?yàn)榇颂幰合嗪蜌庀喟l(fā)生了強(qiáng)烈攪混。隨后，空泡份額迅速上升達(dá)到0.023，而在孔板出口處空泡份額再次突降，這由流通面積突擴(kuò)及出口處氣液兩相的強(qiáng)力攪混兩方面因素造成。之后，由于熱量不斷導(dǎo)入，截面平均空泡份額迅速增加，達(dá)到0.091。

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

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

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

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

4 結(jié)論

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

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

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

參考文獻(xiàn)：

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