耦合一、二次側(cè)換熱的蒸汽發(fā)生器二次側(cè)流場分析

叢騰龍，田文喜，秋穗正，蘇光輝

(西安交通大學(xué) 動(dòng)力工程多相流國家重點(diǎn)實(shí)驗(yàn)室 核科學(xué)與技術(shù)系，陜西 西安 710049)

蒸汽發(fā)生器(SG)作為壓水堆的關(guān)鍵設(shè)備，同時(shí)承擔(dān)著一、二回路側(cè)換熱以及一回路壓力邊界的作用，在反應(yīng)堆運(yùn)行過程中，需保持傳熱管的完整性；然而，在反應(yīng)堆運(yùn)行期間，傳熱管容易出現(xiàn)震動(dòng)、疲勞以及由此帶來的磨損甚至破裂。傳熱管的失效主要是由管束震動(dòng)引起的磨損、管束熱應(yīng)力損壞以及傳熱管污垢沉積腐蝕造成的，對以上因素進(jìn)行分析，需以SG二次側(cè)的速度場、溫度場作為輸入?yún)?shù)。由于SG結(jié)構(gòu)復(fù)雜、體積龐大，對SG二次側(cè)流場進(jìn)行實(shí)驗(yàn)和精細(xì)網(wǎng)格的CFD計(jì)算難以實(shí)現(xiàn)，因此本文采用多孔介質(zhì)模型對SG二次側(cè)流場進(jìn)行數(shù)值分析。

目前已有一些針對SG二次側(cè)流體熱工水力特性的CFD研究，但這些研究仍有待改進(jìn)。在Ferng等[1-3]的研究中，存在如下問題：1) 進(jìn)出口邊界設(shè)置不合理；2) 兩相流計(jì)算關(guān)系式過于粗糙；3) 采用不正確的阻力和湍流計(jì)算模型；4) 計(jì)算忽略下降段、支承板和汽水分離器；5) 二次側(cè)熱源由RELAP5提供，而非耦合計(jì)算。這一系列簡化或不合理假設(shè)，導(dǎo)致計(jì)算結(jié)果與實(shí)際運(yùn)行情況差別很大，計(jì)算結(jié)果在工程應(yīng)用上不具備參考價(jià)值。

基于以上分析，本文采用多孔介質(zhì)方法，對壓水堆SG二次側(cè)流場進(jìn)行數(shù)值模擬，同時(shí)耦合一、二次側(cè)三維換熱。在模擬中，綜合考慮了管束、下降段、汽水分離器、支承板的阻力影響。對于汽水分離器，只考慮其阻力效應(yīng)，不考慮汽水分離。

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

常用的兩相流計(jì)算模型包含均勻流模型、漂移流模型和兩流體模型。均勻流模型將兩相流體等效為單相流體，不能描述兩相間的滑移速度等量。兩流體模型可準(zhǔn)確描述兩相間的熱力和水力學(xué)不平衡性，但由于管束外流動(dòng)的相間作用力以及相間傳熱等模型目前還不完備，采用該模型需引入較多假設(shè)，導(dǎo)致兩流體模型的計(jì)算精度不能保證。漂移流模型將兩相流體視為混合相，但在相間引入漂移速度，考慮相間的水力學(xué)不平衡性，該模型雖從理論上相對于兩流體模型精度較差，但由于該模型需引入的假設(shè)條件較少，計(jì)算精度可滿足SG內(nèi)多孔介質(zhì)模擬的要求[4]?；谝陨峡紤]，本文采用漂移流模型。多孔介質(zhì)內(nèi)的漂移流模型控制方程如下。

質(zhì)量守恒方程：

β·(βρmvm)=0

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

β·(βρmvmvm)=-β

能量守恒方程：

β·βαkvk(ρkEk+p))=

空泡份額方程：

β·(βαgρgvm)=

式中：下標(biāo)m和g分別代表混合物和氣相；β為孔隙率，即網(wǎng)格中流體體積與網(wǎng)格總體積之比；ρ為密度；vm為混合物速度；μm,eff為混合物有效黏度，等于混合物分子黏度和湍流黏度之和；αk為第k相體積份額；p為壓力；S為附加源項(xiàng)，下標(biāo)v、E、m分別代表動(dòng)量、能量、質(zhì)量；keff為有效熱傳導(dǎo)系數(shù)；vdr,k為第k相的漂移速度，漂移速度與相分布計(jì)算關(guān)系式參見文獻(xiàn)[5]。

動(dòng)量、能量、空泡份額方程中的源項(xiàng)定義如下。

對于動(dòng)量方程，其源項(xiàng)表示由于控制體內(nèi)蒸汽發(fā)生器構(gòu)件的存在導(dǎo)致的附加阻力，包括下降段阻力、管束阻力、汽水分離器阻力、支承板阻力。其中，前3個(gè)阻力均為分布阻力，以分布阻力的形式添加到控制體網(wǎng)格，支承板阻力為局部阻力，以多孔階躍模型的方式添加到支承板所在位置的網(wǎng)格界面上。對順流和橫流阻力分別采用管束外順流[6]和橫流[7]阻力經(jīng)驗(yàn)關(guān)系式計(jì)算。支承板、下降段和汽水分離器的阻力，根據(jù)蒸汽發(fā)生器設(shè)計(jì)參數(shù)給定[4]。

對于能量方程，其源項(xiàng)為一次側(cè)向二次側(cè)釋熱。將一次側(cè)流場劃分為與二次側(cè)相同的網(wǎng)格，對于每個(gè)網(wǎng)格，每次迭代開始之前，在UDF中分別計(jì)算網(wǎng)格內(nèi)一次側(cè)向二次側(cè)的釋熱量。U型管內(nèi)為單相液體對流，采用Dittus-Boelter公式[8]計(jì)算；管外換熱包括單相對流換熱、過冷沸騰換熱以及飽和沸騰換熱，對于單相對流換熱，將流速分解為順流流速和橫流流速，分別計(jì)算順流換熱系數(shù)[9-10]和橫流換熱系數(shù)[11]，將這兩者的加權(quán)和[4]作為實(shí)際的斜流換熱系數(shù)。當(dāng)主流溫度達(dá)到飽和溫度時(shí)，發(fā)生飽和沸騰；當(dāng)管外壁溫度大于飽和溫度且主流溫度低于飽和溫度時(shí)，有可能發(fā)生過冷沸騰，本文在該工況下分別按照單相流動(dòng)和過冷沸騰計(jì)算壁面熱流密度，若采用過冷沸騰公式計(jì)算得到的熱流密度大于采用單相流動(dòng)的計(jì)算值，認(rèn)為過冷沸騰發(fā)生。

在計(jì)算中，給定一、二次側(cè)入口的質(zhì)量流速和溫度，即可計(jì)算一、二次側(cè)之間的換熱。

對于空泡份額方程，其源項(xiàng)表示氣相質(zhì)量源項(xiàng)，其值等于蒸發(fā)率與冷凝率之差。動(dòng)量方程中的漂移速度采用Lellouche-Zolotar公式[5]計(jì)算。

對于湍流模型，已公開發(fā)表的模型多是針對純流體區(qū)域的精細(xì)網(wǎng)格提出的[12]，對于多孔介質(zhì)內(nèi)的湍流模型，研究者在多孔介質(zhì)粗網(wǎng)格內(nèi)對已有湍流模型(如k-ε)進(jìn)行積分，得到適用于多孔介質(zhì)的湍流模型[13-15]，但這些模型僅用于單相流動(dòng)，并無針對兩相流動(dòng)的湍流模型。本文采用文獻(xiàn)[16-17]提出的針對管束間兩相流動(dòng)的零方程湍流模型計(jì)算湍流黏度及有效黏度，該模型在ATHOS程序中已被驗(yàn)證[4,18]。

以上模型構(gòu)成完備的控制方程組，將這些方程在ANSYS FLUENT求解器中求解，求解方法為COUPLE算法。

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

本文所選取的二次側(cè)計(jì)算區(qū)域?yàn)椋簭南陆刀稳肟谙蛳碌絿迦笨冢缓笳哿鬟M(jìn)入管束區(qū)，最終進(jìn)入主汽水分離器。計(jì)算區(qū)域示意圖如圖1a所示。對于圓筒形的汽水分離器筒體，為簡化網(wǎng)格劃分，采用方形通道代替，保持筒體流通面積和中心位置不變。為了提高計(jì)算的收斂性并模擬汽水分離器的阻力對上游流場的影響，將汽水分離器添加到計(jì)算域中，汽水分離器只考慮其阻力，不考慮汽水分離效果。對計(jì)算區(qū)域劃分網(wǎng)格，所有網(wǎng)格均為六面體網(wǎng)格，在支承板位置添加多孔階躍邊界模擬支承板阻力。網(wǎng)格劃分如圖1b所示，網(wǎng)格質(zhì)量大于0.5。圖1c為計(jì)算區(qū)域管束布置示意圖。為進(jìn)行網(wǎng)格獨(dú)立性驗(yàn)證，分別對4組網(wǎng)格進(jìn)行計(jì)算，獲得獨(dú)立解時(shí)的網(wǎng)格數(shù)量為111 099，如圖2所示。

a——計(jì)算區(qū)域示意圖；b——網(wǎng)格劃分示意圖；c——管束示意圖

圖2 網(wǎng)格獨(dú)立性驗(yàn)證

二次側(cè)網(wǎng)格入口邊界設(shè)置為速度入口，出口邊界設(shè)置為壓力出口，考慮到多孔介質(zhì)網(wǎng)格尺度較大，壁面無滑移邊界不適用，因此本文采用滑移壁面，將壁面摩擦阻力分布添加到壁面網(wǎng)格內(nèi)。入口溫度由疏水和新給水的流量及焓確定，出口壓力設(shè)置為蒸汽發(fā)生器頂部腔室壓力。

一次側(cè)流場和一、二次側(cè)流體的換熱均采用UDF進(jìn)行計(jì)算，二次側(cè)流場計(jì)算每迭代一步，UDF計(jì)算1次一次側(cè)流場和一、二次側(cè)換熱量，將換熱量作為二次側(cè)流場內(nèi)熱源賦給FLUENT求解器。

3 結(jié)果分析

給定一、二次側(cè)入口流量和溫度，采用FLUENT計(jì)算二次側(cè)流場，同時(shí)采用UDF耦合一次側(cè)流場及一、二次側(cè)換熱，得到三維的一、二次側(cè)流場，管內(nèi)外換熱系數(shù)及U型管溫度分布。FLUENT計(jì)算值與設(shè)計(jì)值的對比列于表1。

表1 FLUENT計(jì)算值與設(shè)計(jì)值的對比

3.1 一、二次側(cè)溫度分析

圖3為SG二次側(cè)流體在對稱截面和水平切面上的溫度云圖，圖3a中右側(cè)(x>0)為熱側(cè)，左側(cè)為冷側(cè)。由圖可看出，熱側(cè)流體在進(jìn)入管束區(qū)后，很快達(dá)到飽和溫度，而冷側(cè)流體緩慢達(dá)到飽和，冷側(cè)流體流過蒸汽發(fā)生器內(nèi)1/3高度時(shí)才完全達(dá)到飽和。達(dá)到飽和后，流體溫度保持不變。由圖3b所示的橫切面溫度分布可看出，由于流動(dòng)和幾何的對稱性，流體溫度在截面上對稱。

a——對稱截面上溫度分布；b——不同高度水平截面上的溫度分布

圖4為徑向積分后溫度沿管束方向的分布。一次側(cè)溫度沿著流動(dòng)的進(jìn)行逐漸降低。二次側(cè)溫度在冷熱側(cè)均隨高度的增大而增大，并最終維持在飽和溫度。熱側(cè)流體上升到約1 m高度時(shí)，即達(dá)到飽和溫度；冷側(cè)流體在上升到約4 m高度時(shí)才達(dá)到飽和溫度。由于二次側(cè)換熱系數(shù)遠(yuǎn)大于一次側(cè)換熱系數(shù)，管外壁溫度沿著管束方向基本不變，維持在547 K。管內(nèi)壁溫度沿著管束方向自熱側(cè)向冷側(cè)逐漸降低，管內(nèi)外壁溫差沿管束方向逐漸減小，管內(nèi)壁最高溫度約為570.6 K。

圖4 溫度沿管束方向的分布

圖5為一、二次側(cè)換熱系數(shù)及壁面熱流密度沿管束的分布。管內(nèi)為單相強(qiáng)制對流，換熱系數(shù)隨物性的變化有輕微下降趨勢；管外分別為單相對流、過冷沸騰和飽和沸騰，其中沸騰換熱系數(shù)遠(yuǎn)大于單相換熱系數(shù)。一次側(cè)平均換熱系數(shù)為15 856.5 W/(m2·K)，二次側(cè)平均換熱系數(shù)為63 623.0 W/(m2·K)，二次側(cè)最大換熱系數(shù)為122 862.9 W/(m2·K)。由于換熱系數(shù)和管內(nèi)外溫差的降低，管壁面熱流密度沿管束方向也逐漸降低，外壁面平均熱流密度為149.9 kW/m2，最大、最小壁面熱流密度分別為369.6和44.6 kW/m2。

圖5 一、二次側(cè)換熱系數(shù)及壁面熱流密度沿管束方向的分布

3.2 二次側(cè)流動(dòng)含汽率分析

SG二次側(cè)流體的流動(dòng)含汽率分布如圖6所示?？煽闯?，在管束區(qū)入口區(qū)域，冷熱側(cè)流動(dòng)含汽率差別較大，熱側(cè)流動(dòng)含汽率小于0.01的區(qū)域很小，而冷側(cè)該區(qū)域占有較大面積。隨著流體沿管束上升，流體不斷蒸發(fā)，氣泡在橫向擴(kuò)散，冷熱側(cè)流動(dòng)含汽率逐漸均衡，且熱側(cè)靠近中心處流動(dòng)含汽率逐漸變?yōu)樽畲蟆?/p>

a——對稱截面上的流動(dòng)含汽率分布；b——不同高度水平截面上的流動(dòng)含汽率分布

圖7為流動(dòng)含汽率沿高度方向的分布，3條曲線分別表示整個(gè)二次側(cè)區(qū)域水平截面平均值、熱側(cè)平均值、冷側(cè)平均值。由圖7可看出，隨高度的增加，流動(dòng)含汽率逐漸增大。熱側(cè)流動(dòng)含汽率遠(yuǎn)大于冷側(cè)流體。冷、熱側(cè)出口平均流動(dòng)含汽率分別為0.16和0.46，總體出口平均流動(dòng)含汽率為0.27。在流體進(jìn)入汽水分離器前，流動(dòng)含汽率出現(xiàn)降低，這是由于在兩相流體流出管束區(qū)、流入汽水分離器前，流通截面積突變，導(dǎo)致氣相流速降低，從而導(dǎo)致流動(dòng)含汽率降低。在流體流入汽水分離器之后，流動(dòng)含汽率保持恒定。

圖7 流動(dòng)含汽率沿高度方向的分布

圖8 汽水分離器內(nèi)流動(dòng)含汽率分布

圖8為汽水分離器內(nèi)流動(dòng)含汽率的分布。由圖8可看出，對于熱側(cè)最內(nèi)層汽水分離器，其流動(dòng)含汽率很高，達(dá)0.75，而對于冷側(cè)外圍的汽水分離器，流動(dòng)含汽率只有0.07，不同的汽水分離器內(nèi)混合物流動(dòng)情況相差很大。

3.3 二次側(cè)壓力分析

圖9為對稱截面上的壓力分布，壓力選取出口截面為參考面，該截面壓力設(shè)為0 Pa，圖中標(biāo)注壓力為相對壓力。在下降段，由于重位壓降、下降腔室壁面摩擦壓降和加速壓降的共同作用，壓力沿著下降段向下逐漸增大；流體在圍板缺口區(qū)域，由于流體轉(zhuǎn)向，壓力陡增，之后進(jìn)入管束區(qū)域，流體壓力隨著高度方向逐漸減小。同一高度截面上，熱側(cè)壓力總體稍大于冷側(cè)壓力。圖10為壓力沿高度方向的分布。壓力沿著高度方向出現(xiàn)明顯的階躍變化，這是由于支承板的存在，引入局部阻力，導(dǎo)致流體經(jīng)過支承板時(shí)，壓力突降；支承板導(dǎo)致壓降的大小隨軸向高度的增大而增大，這是由于隨高度的增加，兩相流體的流動(dòng)含汽率增大，兩相摩擦壓降倍增因子增大，從而在相同的阻力系數(shù)下，支承板位置越高，引入的壓降越大。位于11～12 m之間的壓力階躍是由于流體流入汽水分離器造成的，流體流通截面積突降，流速增大，導(dǎo)致靜壓降低。

圖9 壓力在對稱截面上的分布

圖10 壓力沿高度方向的分布

3.4 二次側(cè)速度分析

圖11為對稱截面上的速度分布。由圖11可看出，熱側(cè)流體速度明顯大于冷側(cè)，最大速度約7.4 m/s。還可看出，除圍板缺口折流部分和彎管區(qū)域，流體基本沿管束方向流動(dòng)，橫流較小。

圖11 對稱截面上的速度矢量圖

圖12 彎管部分橫流速度分布及橫流能量分布

圖12為彎管部分橫流速度分布及橫流能量分布。由圖12可看出，自冷側(cè)到熱側(cè)，橫流速度整體呈現(xiàn)先增大后減小的趨勢，最大橫流速度約為4.06 m/s；在彎管部分頂端，由于汽水分離器的影響，導(dǎo)致汽水分離器控制部分的流速高于彎管最頂端的速度，這與Ferng等[2-3]的不考慮汽水分離器的速度分布不同。根據(jù)文獻(xiàn)[20]，管束流致振動(dòng)強(qiáng)度與橫流能量(ρu2)呈正比，由于冷側(cè)流體密度較大，冷側(cè)橫流能量高于熱側(cè)，冷熱側(cè)最大橫流能量分別為1 145和695 J/m3。

4 結(jié)語

本文采用多孔介質(zhì)模型，對壓水堆蒸汽發(fā)生器二次側(cè)三維兩相流場進(jìn)行數(shù)值模擬，同時(shí)耦合一、二次側(cè)三維換熱，得到如下結(jié)論：

1) 蒸汽發(fā)生器內(nèi)冷熱側(cè)流體由于受熱不均，導(dǎo)致兩側(cè)流體流動(dòng)狀態(tài)差別較大；

2) 一、二次側(cè)平均換熱系數(shù)分別為15 856.5和63 623.0 W/(m2·K)，二次側(cè)最大換熱系數(shù)為122 862.9 W/(m2·K)，U型管外壁面平均熱流密度為149.9 kW/m2；

3) 汽水分離器內(nèi)流動(dòng)含汽率分布極不均勻，其值介于0.07～0.75之間，此分布可為汽水分離器設(shè)計(jì)提供入口邊界條件；

4) U型管彎管段最大橫流速度約為4.06 m/s；冷側(cè)沖刷U型管的橫流能量大于熱側(cè)，其最大值分別為1 145和695 J/m3。汽水分離器的位置對U型管橫流速度以及橫流能量的影響很大。

參考文獻(xiàn)：

[1] FERNG Y M. Investigating the distribution characteristics of boiling flow and released nuclide in the steam generator secondary side using CFD methodology[J]. Ann Nucl Energy, 2007, 34(9): 724-731.

[2] FERNG Y M, CHANG H J. CFD investigating the impacts of changing operating conditions on the thermal-hydraulic characteristics in a steam generator[J]. Appl Therm Eng, 2008, 28(5): 414-422.

[3] FERNG Y M, YINPANG M, KANG J C. Thermal-hydraulic simulation of localized flow characteristics in a steam generator[J]. Nucl Technol, 2001, 136(2): 186-196.

[4] KEETON L, SINGHAL A, SRIKANTIAH G. ATHOS3: A computer program for thermal-hydraulic analysis of steam generators[M]. Palo Alto, CA, US: Electric Power Research Institute, 1986.

[5] LELLOUCHE G, ZOLOTAR B. A mechanistic model for predicting two-phase void for water in vertical tubes, channels and rod bundles, EPRI-N-2246[R]. Palo Alto, CA, US: Electric Power Research Institute, 1982.

[6] MacADAMS W H. Heat transmission[M]. New York, US: McGraw Hill, 1954.

[7] GRIMISON E. Correlation and utilization of new data on flow resistance and heat transfer for cross flow of gases over tube banks[J]. Trans ASME, 1937, 59(7): 583-594.

[8] DITTUS F W, BOELTER L M K. Heat transfer in automobile radiators of the tubular type[J]. Int Commun Heat Mass Transfer, 1985, 12(1): 3-22.

[9] DINGEE D A, BELL W, CHASTAIN J, et al. Heat transfer from parallel rods in axial flow[R]. Columbus, Ohio: Battelle Memorial Inst., 1955.

[10] DINGEE D A, CHASTAIN J W. Heat transfer from parallel rods in an axial flow[C]∥ASME Reactor Heat Transfer Conference. Oak Ridge: US Atomic Energy Commission, 1968.

[11] DWYER O, SHEEHAN T, WEISMAN J, et al. Cross flow of water through a tube bank at reynolds numbers up to a million[J]. Industrial & Engineering Chemistry, 1956, 48(10): 1 836-1 846.

[12] 陶文銓. 數(shù)值傳熱學(xué)[M]. 第2版. 西安：西安交通大學(xué)出版社，2001.

[13] ANTOHE B, LAGE J. A general two-equation macroscopic turbulence model for incompressible flow in porous media[J]. Int J Heat Mass Transfer, 1997, 40(13): 3 013-3 024.

[14] CHANDESRIS M, SERRE G, SAGAUT P. A macroscopic turbulence model for flow in porous media suited for channel, pipe and rod bundle flows[J]. Int J Heat Mass Transfer, 2006, 49(15): 2 739-2 750.

[15] De LEMOS M J S. Turbulence in porous media: Modeling and applications[M]. London: Elsevier, 2012.

[16] LAUNDER B E, SPALDING D B. Lectures in mathematical models of turbulence[M]. London: Academic Press, 1972.

[17] SCHLICHTING H. Boundary-layer theory[M]. 7th ed. Germany: Springer, 1979.

[18] HOPKINS G. Verification of the ATHOS3 code against feedring and preheat steam generator test data[M]. Palo Alto, CA, US: Electric Power Research Institute, 1988.

[19] 林誠格，郁祖盛，歐陽予. 非能動(dòng)安全先進(jìn)核電廠AP1000[M]. 北京：原子能出版社，2008.

[20] AXISA F, ANTUNES J, VILLARD B. Overview of numerical methods for predicting flow-induced vibration[J]. J Pressure Vessel Technol, 1988, 110(1): 6-14.