蒸汽發(fā)生器二回路過(guò)冷沸騰的數(shù)值模擬*

張小英 丁斐 陳佳躍

(華南理工大學(xué) 電力學(xué)院，廣東 廣州 510640)

蒸汽發(fā)生器是壓水堆核電站中連接一、二回路的關(guān)鍵設(shè)備，也是核電廠工作運(yùn)行中故障發(fā)生最多的設(shè)備之一.在蒸汽發(fā)生器二回路的預(yù)熱段易出現(xiàn)過(guò)冷沸騰.由于過(guò)冷沸騰區(qū)域內(nèi)不穩(wěn)定的兩相流動(dòng)徑向存在較大的溫差，導(dǎo)致兩相之間出現(xiàn)熱力學(xué)與動(dòng)力學(xué)不平衡，有可能引發(fā)自然循環(huán)系統(tǒng)不穩(wěn)定性.文獻(xiàn)［1］通過(guò)實(shí)驗(yàn)研究了自然循環(huán)系統(tǒng)中過(guò)冷沸騰流動(dòng)的不穩(wěn)定性，文獻(xiàn)［2］基于RELAP5 程序分析了單通道內(nèi)自然循環(huán)流動(dòng)的不穩(wěn)定性.然而，當(dāng)前已公布的關(guān)于蒸汽發(fā)生器二次側(cè)流體的熱工水力分析一般不考慮預(yù)熱段的過(guò)冷沸騰［3-6］，或簡(jiǎn)化為單相換熱處理［7］.為優(yōu)化蒸汽發(fā)生器運(yùn)行穩(wěn)定性，對(duì)蒸汽發(fā)生器二回路的過(guò)冷沸騰的研究是很有必要的.

Saha 模型和Levy 模型是目前國(guó)際上應(yīng)用最為廣泛的過(guò)冷沸騰模型.但孫奇等［8-9］通過(guò)實(shí)驗(yàn)證明，在低流速過(guò)冷沸騰下，Saha 模型和Levy 模型預(yù)測(cè)結(jié)果均明顯高于實(shí)驗(yàn)值，因此兩種模型不能準(zhǔn)確預(yù)測(cè)低流速下的過(guò)冷沸騰.在蒸汽發(fā)生器二回路預(yù)熱段流體的流速一般較低，不能采用現(xiàn)有的過(guò)冷沸騰計(jì)算模型對(duì)其進(jìn)行模擬.RPI 模型通過(guò)計(jì)算氣泡的參量，基于熱量平衡而求得壁面氣化量，被廣泛地應(yīng)用于過(guò)冷沸騰區(qū)內(nèi)對(duì)壁面熱量的劃分［10-12］.

目前CFD 程序已經(jīng)被廣泛地用于兩相領(lǐng)域的研究［12-14］.為了充分認(rèn)識(shí)蒸汽發(fā)生器二次側(cè)預(yù)熱段低流速情況下過(guò)冷沸騰現(xiàn)象和機(jī)理，文中基于CFD 程序現(xiàn)有的兩流體模型，嵌入RPI 過(guò)冷沸騰模型對(duì)其修正和簡(jiǎn)化，并將該模型的模擬結(jié)果與Bartolomej 等［15］的實(shí)驗(yàn)數(shù)據(jù)對(duì)比，驗(yàn)證了RPI 模型用于低流速過(guò)冷沸騰模擬問(wèn)題的適應(yīng)性;采用修正后的兩流體模型對(duì)大亞灣900 MW 壓水堆核電站蒸汽發(fā)生器二回路的過(guò)冷沸騰進(jìn)行預(yù)測(cè)，得出蒸汽發(fā)生器二回路管束通道內(nèi)過(guò)冷沸騰的流動(dòng)參數(shù)特征.

1 理論模型

1.1 相間熱量與質(zhì)量的交換

在過(guò)冷沸騰起始點(diǎn)，主流溫度仍低于流體的飽和溫度，流道壁面開(kāi)始出現(xiàn)氣泡時(shí)，液相向氣相轉(zhuǎn)化;隨著流體繼續(xù)被加熱至充分發(fā)展段起始點(diǎn)后，氣泡開(kāi)始脫離壁面，進(jìn)入欠熱主流流體從而被冷凝，氣相又向液相轉(zhuǎn)化.因此，在過(guò)冷沸騰區(qū)內(nèi)，液相轉(zhuǎn)化為氣相的同時(shí)也存在著氣相向液相的轉(zhuǎn)化，兩相之間熱量與質(zhì)量的傳遞包括由對(duì)流引起的傳遞和相變引起的傳遞.

1.1.1 質(zhì)量的交換

在過(guò)冷沸騰換熱中，兩相間的質(zhì)量傳遞包括液相的蒸發(fā)和氣相的冷凝所引起的變化.液相的蒸發(fā)量與生成氣泡的參數(shù)有關(guān)，氣相的冷凝量主要取決于流體的液相溫度，文中按照文獻(xiàn)［11］中推薦的公式計(jì)算，具體如下.

(1)液相的蒸發(fā)量

式中，NW為汽化核心密度，f 為氣泡脫離頻率，ρg為氣相密度，VW為氣泡脫離體積.其中:

DW采用Tolubinsky 關(guān)系式得到:

式中:TW為壁面溫度;Tl為液相溫度;g 為重力加速度;ρl為液相密度;DW為氣泡脫離直徑;ΔTsub為壁面過(guò)熱度.

(2)氣相的冷凝量

式中，Tsat為飽和溫度，Ai為受熱面積，Hlg為氣化潛熱，hlg為換熱系數(shù)，hlg可按照Ranz_Marshll 關(guān)聯(lián)式計(jì)算，如下:

式中，l為流體液相的導(dǎo)熱系數(shù)，Re 為雷諾數(shù)，Pr 為普朗特?cái)?shù)，α 為空泡份額.

1.1.2 動(dòng)量的交換

過(guò)冷沸騰中動(dòng)量的交換主要考慮由于界面作用力引起的傳遞.界面作用力為

式中:FD為拖曳力;Flift為升力;FTD為湍流耗散力.各項(xiàng)界面作用力按照文獻(xiàn)［9］中的推薦公式計(jì)算，如下:

式中，F(xiàn)TD，l、FTD，g分別為液、氣相的湍流耗散力，CD、Cl、CTD分別為曳力系數(shù)、升力系數(shù)和湍流耗散系數(shù)，vg、vl分別為氣液兩相的速度，ωg為氣相速度旋度，k為湍流脈動(dòng)動(dòng)能，為梯度算子.

1.1.3 能量的交換

過(guò)冷沸騰中兩相之間能量的傳遞可包括由于相變引起的能量交換和由于兩相溫差引起的對(duì)流換熱.相變引起的能量交換［11］為

兩相溫差引起的對(duì)流換熱可按文獻(xiàn)［7］中推薦的公式計(jì)算:

式中相間換熱系數(shù)hgl、面積Ai可按式(6)、(7)計(jì)算.

1.2 RPI 過(guò)冷沸騰模型

在過(guò)冷沸騰區(qū)域內(nèi)，壁面?zhèn)鳠岬挠?jì)算十分重要，文中選擇Kurul 等［16］提出的RPI 模型，該模型基于熱量平衡，把過(guò)冷沸騰區(qū)壁面與流體之間的換熱量QW分為三部分——液相對(duì)流換熱帶走的熱量QC、氣泡脫離壁面后中心過(guò)冷流體填充其位置使壁面驟冷而交換的熱量QQ和液相蒸發(fā)而帶走的熱量QE:

(1)對(duì)流換熱量為

式中:St 為局部斯坦頓數(shù)，St =Nu/(Re×Pr);cp，l為液相的比定壓熱容;氣相面積Ab=min(1，NWDW2/4);汽化核心密度NW、氣泡脫離直徑DW可按式(2)、(4)計(jì)算.

(2)驟冷換熱量為

式中，t 為氣泡脫離周期，t=0.8/f，f 可由式(3)計(jì)算得到.

(3)蒸發(fā)換熱量為

其中氣泡的體積VW可根據(jù)氣泡脫離直徑DW計(jì)算.

2 算例計(jì)算

2.1 模型驗(yàn)證算例

為驗(yàn)證RPI 模型的適應(yīng)性，文中對(duì)Bartolomej等［15］研究的豎直圓管高壓過(guò)冷沸騰實(shí)驗(yàn)的熱工參數(shù)進(jìn)行了計(jì)算，并將計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比.計(jì)算條件為:管徑15.4mm，加熱段長(zhǎng)度2m，壁面熱流量5.7×105W/m2，入口質(zhì)量流速900kg/(m2·s)(等效為1.143 m/s，蒸汽發(fā)生器預(yù)熱段流速為1 m/s)，入口過(guò)冷度為64.5 K，工作壓力為4.5 MPa.

采用RPI 模型計(jì)算所得的空泡份額(α)、采用Levy 模型計(jì)算的結(jié)果與Bartolomej 實(shí)驗(yàn)結(jié)果見(jiàn)圖1，其中Z 為豎直圓管的豎直方向.可以看出采用RPI模型計(jì)算所得的空泡份額與實(shí)驗(yàn)值相符較好，而Levy 模型所得的空泡份額明顯高于實(shí)驗(yàn)值，可見(jiàn)RPI 模型更加適用于對(duì)高壓、低流速下過(guò)冷沸騰換熱特性的預(yù)測(cè).

圖1 豎直圓管內(nèi)平均空泡份額的分布Fig.1 Distributions of average void fraction in vertical tube

2.2 實(shí)例計(jì)算

文中以大亞灣900 MW 核電站蒸汽發(fā)生器為例，選取蒸汽發(fā)生器管束間的一個(gè)垂直流道單元作為研究對(duì)象，見(jiàn)圖2.換熱管外徑為19.05 mm，壁厚為1.09 mm，材料為奧氏不銹鋼，節(jié)距為25 mm，流道單元長(zhǎng)1 000 mm.管壁材料主要物性參數(shù)為:密度8 030 kg/m3，導(dǎo)熱系數(shù)18.3 W/(m·K)，比等壓熱容430 J/(kg·K).

圖2 流道示意圖Fig.2 Schematic diagram of flow channel

計(jì)算所選取的邊界條件如下:入口為速度入口邊界，速度為1 m/s，過(guò)冷度為15 K;出口為壓力出口邊界，壓力為0 Pa.如圖2 所示，流道單元的圓弧面為受熱面，給定對(duì)流換熱系數(shù)為41 080 W/(m2·K)(根據(jù)Dittus-Boelter 公式計(jì)算所得)，給定一回路平均溫度583.15 K;其他壁面邊界為對(duì)稱(chēng)邊界，視為絕熱;工作壓力設(shè)為6.89 MPa.假定流體為不可壓縮，流體的物性參數(shù)考慮隨溫度而變化.

通過(guò)用戶(hù)自定義函數(shù)，采用RPI 模型對(duì)CFD 程序中現(xiàn)有兩流體模型進(jìn)行了修正，以模擬高壓、低流速下的過(guò)冷沸騰換熱.在求解計(jì)算中，采用有限容積法對(duì)方程組進(jìn)行離散，湍流模型選用RNG k-ε 模型，對(duì)壓力和速度的耦合計(jì)算選用SIMPLE 算法，當(dāng)計(jì)算殘差、進(jìn)出口相對(duì)質(zhì)量差小于10-4時(shí)可以認(rèn)為計(jì)算收斂.

3 計(jì)算結(jié)果與分析

采用RPI 模型計(jì)算流道中空泡份額的結(jié)果見(jiàn)圖3、4.從圖3 看到:空泡份額開(kāi)始很小，隨后近似呈線性增加，且增加速率逐漸變大.圖4 示出了流體在不同截面的空泡份額分布，可以看出壁面附近空泡份額明顯高于主流體區(qū)域.原因是初始段流體的過(guò)冷度較大，流體在壁面附近先生成氣泡，隨著氣泡的脫離和破裂，氣體開(kāi)始進(jìn)入主流體，隨著過(guò)冷度的減小，流體溫度逐漸接近飽和溫度，氣化速率逐漸變大.

圖3 平均空泡份額沿軸向變化曲線Fig.3 Curve of the average void fraction along the axial

圖4 不同截面空泡份額分布情況Fig.4 Distributions of void fraction in different cross-sections

圖5 示出了流體的平均溫度(T)在不同截面的分布情況，可以看出:隨著壁面對(duì)流體的加熱，流體溫度逐漸上升，管壁附近流體溫度明顯高于主流的溫度.

圖5 不同截面溫度分布情況Fig.5 Distributions of temperature in different cross-sections

圖6 示出了流體、壁面平均溫度沿軸向的變化，可以看出:流體的平均溫度明顯低于壁面溫度，隨著軸線距離增加，流體與壁面的溫差逐漸減小;在Z =700 mm 時(shí)壁面溫度開(kāi)始超過(guò)飽和溫度，但主流的平均溫度仍低于飽和溫度，表明該區(qū)域出現(xiàn)過(guò)冷沸騰.

圖7 示出了文中采用RPI 模型計(jì)算壁面對(duì)流體加熱熱流的分配情況.可以看到:①在過(guò)冷沸騰初始階段，單相對(duì)流換熱量較大，而蒸發(fā)換熱量和驟冷換熱量較小，這是由于初始段氣相比例很小，熱量的傳遞重要是單相的對(duì)流換熱;②隨著主流溫度的升高，空泡份額逐漸增加，導(dǎo)致對(duì)流換熱量下降，蒸發(fā)換熱量上升;③氣泡脫離壁面時(shí)，主流流體通過(guò)填充其位置而使壁面驟冷，隨著氣泡的大量生成，壁面附近的液相流體越來(lái)越少，因此驟冷換熱量呈現(xiàn)一個(gè)先上升而后平緩下降的過(guò)程.

圖6 溫度沿軸向變化曲線Fig.6 Curve of the temperature along the axial

圖7 壁面熱流量分配情況Fig.7 Distributions of heat flux in wall

圖8、9 分別示出了流體在不同截面液相、氣相速度的分布情況.可以看出:①在流道的下游，氣液兩相速度均增大，超過(guò)入口速度，這是由于隨著氣相的生成，流體的平均密度降低，而流道內(nèi)流體的體積流量不變，導(dǎo)致流體的整體速度增大;②液相速度在通道入口端變化不大，在出口端隨著空泡份額增加，壓力逐漸下降，液相速度逐漸變大;氣相速度明顯高于液相的速度，是由于氣相不僅在氣液界面剪切力的作用下隨液相一起運(yùn)動(dòng)，還受到浮升力的作用，氣相溫度越高，空泡份額越大，所受到的浮力越大，速度也就越大;③在Z =250 mm 處，壁面附近區(qū)域的液相流速很低，是由于壁面剛開(kāi)始產(chǎn)生氣泡，氣泡尚未脫離壁面，阻礙了流體的流動(dòng).

圖8 不同截面液相速度分布情況Fig.8 Distributions of liquid speed in different cross-sections

圖9 不同截面氣相速度分布情況Fig.9 Distributions of vapor speed in different cross-sections

4 結(jié)論

文中采用RPI 模型對(duì)CFD 程序現(xiàn)有的兩流體模型進(jìn)行修正，并與Bartolomej 等的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，發(fā)現(xiàn):與Levy 模型相比，采用修正后的模型計(jì)算預(yù)測(cè)段過(guò)冷沸騰所得的空泡份額與實(shí)驗(yàn)值符合得較好.

以大亞灣核電站為例，采用RPI 模型和CFD 程序模擬蒸汽發(fā)生器二回路預(yù)熱段的過(guò)冷沸騰，得出了單元通道內(nèi)流體的空泡份額、平均溫度沿著通道方向的分布.結(jié)果表明:壁面熱量的分布沿著流動(dòng)方向有明顯變化，沿流道管長(zhǎng)對(duì)流換熱量逐漸下降，蒸發(fā)換熱量迅速上升，驟冷換熱量則呈現(xiàn)先上升而后平緩下降的趨勢(shì).本研究方法和結(jié)果對(duì)蒸汽發(fā)生器二回路預(yù)熱段的沸騰換熱研究具有一定的參考意義，進(jìn)一步的研究還需要考慮管束間多個(gè)通道的影響情況.

［1］楊瑞昌，唐虹，王彥武.自然循環(huán)過(guò)冷沸騰流動(dòng)不穩(wěn)定性的研究［J］.工程熱物理學(xué)報(bào)，2004，25(3):435-438.Yang Rui-Chang，Tang Hong，Wang Yan-wu.Study on flow instability in a natural circulation system with subcooled boiling ［J］.Journal of Engineering Thermophysics，2004，25(3):435-438.

［2］刑立淼，郭赟，曾和義.基于RELAP5 的單通道自然循環(huán)流動(dòng)不穩(wěn)定性分析［J］.原子能科學(xué)技術(shù)，2010，44(8):958-963.Xing Li-miao，Guo Yun，Zeng He-yi.Investigation on circulation flow instability in single channel based on RELAP5 code［J］.Atomic Energy Science and Technology，2010，44(8):958-963.

［3］楊晨，唐勝利，何祖威.U 型管蒸汽發(fā)生器數(shù)學(xué)模型的建立及其仿真研究［J］.系統(tǒng)仿真學(xué)報(bào)，1999，11(3):168-170，193.Yang Chen，Tang Sheng-li，He Zu-wei.Development of Utube steam generator model and its simulation［J］.Journal of System Simulation，1999，11(3):168-170，193.

［4］Ferng Y M，Ma Y P，Kang J C.Thermal-hydraulic simulation of localized flow characteristics in a steam generator［J］Nuclear Technology，2001，136(2):186-196.

［5］Ferng Y M，Chang H J.CFD investigating the impacts ofchanging operating conditions on the thermal-hydraulic characteristics in a steam generator［J］.Applied Thermal Engineering，2008，28(5/6):414-422.

［6］Ferng Y M.Investigating the distribution characteristics of boiling flow and released nuclide in the steam generator secondary side using CFD methodology ［J］.Annals of Nuclear Energy，2007，34(9):724-731.

［7］張龍飛，張大發(fā).U 型管蒸汽發(fā)生器數(shù)學(xué)模型及動(dòng)態(tài)仿真［J］.鍋爐技術(shù)，2006，37(4):9-11，39.Zhang Long-fei，Zhang Da-fa.U-tube steam generator model and dynamic simulation ［J］.Boiler Technology，2006，37(4):9-11，39.

［8］Sun Q，Yang R，Zhao H.Predictive study of the incipient point of net vapor generation in low-flow subcooled boiling［J］.Nuclear Engineering and Design，2003，225(2/3):249-256.

［9］孫奇，楊瑞昌，王小軍，等.低流速過(guò)冷沸騰截面平均空泡率實(shí)驗(yàn)研究［J］.原子能科學(xué)技術(shù)，2005，39(3):193-197.Sun Qi，Yang Rui-chang，Wang Xiao-jun，et al.Experimental study of average void fraction in low-flow subcooled boiling［J］.Atomic Energy Science and Technology，2005，39(3):193-197.

［10］Chen E，Li Y，Cheng X.CFD simulation of upward subcooled boiling flow of refrigerant-113 using the two-fluid model［J］.Applied Thermal Engineering，2009，29(11/12):2508-2517.

［11］陳二鋒，厲彥忠，王斯民.豎直環(huán)管內(nèi)低壓水過(guò)冷沸騰數(shù)值模擬［J］.西安交通大學(xué)學(xué)報(bào)，2008，42(7):855-859.Chen Er-feng，Li Yan-zhong，Wang Si-min.Numerical simulation of subcooled boiling water in vertical concentric annulus under low pressure［J］.Journal of Xi’an Jiaotong University，2008，42(7):855-859.

［12］樊普，賈斗南，秋穗正.低壓下水欠熱流動(dòng)沸騰的兩相CFD 數(shù)值模擬研究［J］.原子能科學(xué)技術(shù)，2011，45(4):412-420.Fan Pu，Jia Dou-nan，Qiu Sui-zheng.CFD investigation of subcooled flow boiling model under low pressure［J］.Atomic Energy Science and Technology，2011，45(4):412-420.

［13］朱冬生，張景衛(wèi)，吳治將，等.板式蒸發(fā)式冷凝器兩相降膜流動(dòng)CFD 模擬及傳熱研究［J］.華南理工大學(xué)學(xué)報(bào):自然科學(xué)版，2008，36(7):6-10.Zhu Dong-sheng，Zhang Jing-wei，Wu Zhi-jiang，et al.CFD simulation and investigation into heat transfer for falling film with two-phase flow in plate-type evaporative condenser ［J］.Journal of South China University of Technology:Natural Science Edition，2008，36(7):6-10.

［14］王繼紅，張騰飛，王樹(shù)剛，等.水平管道內(nèi)固液兩相流流動(dòng)特性的CFD 模擬［J］.化工學(xué)報(bào)，2011，62(12):3399-3404.Wang Ji-hong，Zhang Teng-fei，Wang Shu-gang，et al.Numerical simulation of liquid-solid slurry flow in horizontal pipeline ［J］.CIESC Journal，2011，62(12):3399-3404.

［15］Bartolomej G G，Chanturiya V M.Experimental study of true void fraction when boiling subcooled water in vertical tubes［J］.Thermal Engineering，1967，14(2):123-128.

［16］Kurul N，Podowski M Z.On the modeling of multidimensional effects in boiling channels ［C］∥ANS Proceedings of 27th National Heat Transfer Conference.Minneapolis:ANS，1991.