蒸汽發生器二回路過冷沸騰的數值模擬*

張小英 丁斐 陳佳躍

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

蒸汽發生器是壓水堆核電站中連接一、二回路的關鍵設備，也是核電廠工作運行中故障發生最多的設備之一.在蒸汽發生器二回路的預熱段易出現過冷沸騰.由于過冷沸騰區域內不穩定的兩相流動徑向存在較大的溫差，導致兩相之間出現熱力學與動力學不平衡，有可能引發自然循環系統不穩定性.文獻［1］通過實驗研究了自然循環系統中過冷沸騰流動的不穩定性，文獻［2］基于RELAP5 程序分析了單通道內自然循環流動的不穩定性.然而，當前已公布的關于蒸汽發生器二次側流體的熱工水力分析一般不考慮預熱段的過冷沸騰［3-6］，或簡化為單相換熱處理［7］.為優化蒸汽發生器運行穩定性，對蒸汽發生器二回路的過冷沸騰的研究是很有必要的.

Saha 模型和Levy 模型是目前國際上應用最為廣泛的過冷沸騰模型.但孫奇等［8-9］通過實驗證明，在低流速過冷沸騰下，Saha 模型和Levy 模型預測結果均明顯高于實驗值，因此兩種模型不能準確預測低流速下的過冷沸騰.在蒸汽發生器二回路預熱段流體的流速一般較低，不能采用現有的過冷沸騰計算模型對其進行模擬.RPI 模型通過計算氣泡的參量，基于熱量平衡而求得壁面氣化量，被廣泛地應用于過冷沸騰區內對壁面熱量的劃分［10-12］.

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

1 理論模型

1.1 相間熱量與質量的交換

在過冷沸騰起始點，主流溫度仍低于流體的飽和溫度，流道壁面開始出現氣泡時，液相向氣相轉化;隨著流體繼續被加熱至充分發展段起始點后，氣泡開始脫離壁面，進入欠熱主流流體從而被冷凝，氣相又向液相轉化.因此，在過冷沸騰區內，液相轉化為氣相的同時也存在著氣相向液相的轉化，兩相之間熱量與質量的傳遞包括由對流引起的傳遞和相變引起的傳遞.

1.1.1 質量的交換

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

(1)液相的蒸發量

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

DW采用Tolubinsky 關系式得到:

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

(2)氣相的冷凝量

式中，Tsat為飽和溫度，Ai為受熱面積，Hlg為氣化潛熱，hlg為換熱系數，hlg可按照Ranz_Marshll 關聯式計算，如下:

式中，l為流體液相的導熱系數，Re 為雷諾數，Pr 為普朗特數，α 為空泡份額.

1.1.2 動量的交換

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

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

式中，FTD，l、FTD，g分別為液、氣相的湍流耗散力，CD、Cl、CTD分別為曳力系數、升力系數和湍流耗散系數，vg、vl分別為氣液兩相的速度，ωg為氣相速度旋度，k為湍流脈動動能，為梯度算子.

1.1.3 能量的交換

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

兩相溫差引起的對流換熱可按文獻［7］中推薦的公式計算:

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

1.2 RPI 過冷沸騰模型

在過冷沸騰區域內，壁面傳熱的計算十分重要，文中選擇Kurul 等［16］提出的RPI 模型，該模型基于熱量平衡，把過冷沸騰區壁面與流體之間的換熱量QW分為三部分——液相對流換熱帶走的熱量QC、氣泡脫離壁面后中心過冷流體填充其位置使壁面驟冷而交換的熱量QQ和液相蒸發而帶走的熱量QE:

(1)對流換熱量為

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

(2)驟冷換熱量為

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

(3)蒸發換熱量為

其中氣泡的體積VW可根據氣泡脫離直徑DW計算.

2 算例計算

2.1 模型驗證算例

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

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

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

2.2 實例計算

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

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

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

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

3 計算結果與分析

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

圖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)在不同截面的分布情況，可以看出:隨著壁面對流體的加熱，流體溫度逐漸上升，管壁附近流體溫度明顯高于主流的溫度.

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

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

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

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

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

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

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

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

4 結論

文中采用RPI 模型對CFD 程序現有的兩流體模型進行修正，并與Bartolomej 等的實驗數據進行對比，發現:與Levy 模型相比，采用修正后的模型計算預測段過冷沸騰所得的空泡份額與實驗值符合得較好.

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

［1］楊瑞昌，唐虹，王彥武.自然循環過冷沸騰流動不穩定性的研究［J］.工程熱物理學報，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 的單通道自然循環流動不穩定性分析［J］.原子能科學技術，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 型管蒸汽發生器數學模型的建立及其仿真研究［J］.系統仿真學報，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］張龍飛，張大發.U 型管蒸汽發生器數學模型及動態仿真［J］.鍋爐技術，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］孫奇，楊瑞昌，王小軍，等.低流速過冷沸騰截面平均空泡率實驗研究［J］.原子能科學技術，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］陳二鋒，厲彥忠，王斯民.豎直環管內低壓水過冷沸騰數值模擬［J］.西安交通大學學報，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］樊普，賈斗南，秋穗正.低壓下水欠熱流動沸騰的兩相CFD 數值模擬研究［J］.原子能科學技術，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］朱冬生，張景衛，吳治將，等.板式蒸發式冷凝器兩相降膜流動CFD 模擬及傳熱研究［J］.華南理工大學學報:自然科學版，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］王繼紅，張騰飛，王樹剛，等.水平管道內固液兩相流流動特性的CFD 模擬［J］.化工學報，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.