壓水堆蒸汽發生器一、二次側穩態流場耦合分析

叢騰龍，田文喜，秋穗正，蘇光輝，謝永誠，姚彥貴

(1.西安交通大學 能源與動力工程學院，陜西 西安 710049；2.上海核工程研究設計院，上海 200233)

蒸汽發生器(SG)作為壓水堆的關鍵設備同時承擔著一、二回路側換熱以及一回路壓力邊界的作用，在反應堆運行過程中需保持傳熱管的完整性；然而在反應堆運行期間，傳熱管易出現震動、疲勞及由此帶來的磨損甚至破裂。傳熱管的失效主要是由于管束震動引起的磨損、管束熱應力損壞及事故工況時管束最高溫度超限造成的，對以上因素進行分析需以SG二次側的速度場、溫度場作為輸入參數。由于SG結構復雜、體積龐大，對SG二次側流場進行實驗和精細網格的CFD計算難以實現，因此本文采用多孔介質模型對SG二次側流場進行數值分析。

最早將多孔介質模型應用于換熱器的是英國學者Patankar等[1]，他們于20世紀70年代就采用此模型對換熱器殼側進行了數值模擬。之后，Sha等[2]又在文獻[1]的基礎上模擬了SG和反應堆堆芯中冷卻劑的流動，Prithiviraj等[3-4]模擬了三維換熱器中的流動。Ferng等[5-7]采用CFX軟件多孔介質模型對SG管束區的二次側流場進行了分析。

目前雖已有一些針對SG二次側流體熱工水力特性的CFD研究，但這些研究仍有待改進。在Ferng等[5-12]的研究中，存在如下問題：1) 進出口邊界設置不合理；2) 兩相流計算關系式過于粗糙；3) 采用不正確的阻力和湍流計算模型；4) 計算忽略下降段、支承板和汽水分離器，這一系列簡化或不合理假設，導致計算結果和實際運行情況差別很大，計算結果在工程應用上不具備參考價值。基于以上分析，本文采用多孔介質模型，對壓水堆SG二次側流場進行數值模擬，同時耦合一、二次側一維換熱。在模擬中綜合考慮管束、下降段、汽水分離器、支承板的阻力影響。

1 數學物理模型

常用的兩相流計算模型包含均勻流模型、漂移流模型和兩流體模型。均勻流模型將兩相流體等效為單相流體，不能描述兩相間的滑移速度等。兩流體模型可準確描述兩相間的熱力和水力學不平衡性，但由于管束外流動的相間作用力及相間傳熱等模型目前尚不完備，采用該模型需引入較多假設，導致兩流體模型的計算精度不能保證。漂移流模型將兩相流體視為混合相，但在相間引入漂移速度，考慮相間的水力學不平衡性，該模型雖從理論上相對于兩流體模型精度較差，但由于該模型需引入的假設條件較少，計算精度可滿足SG內多孔介質模擬的要求[13]。基于上述考慮本文采用漂移流模型。多孔介質內的漂移流模型控制方程如下。

質量守恒方程：

β·(βρmvm)=0

(1)

β·(βρmvmvm)=-β

(2)

能量守恒方程：

β·

(3)

空泡份額方程：

β·(βαgρgvm)=

(4)

式中：β為孔隙率，即網格中流體體積與網格總體積之比；ρ為密度，其中下標m和g分別代表混合物和氣相；vm為混合物速度；μm,eff為混合物有效黏度，等于混合物分子黏度和湍流黏度之和；αk為第k相體積份額；p為壓力；S為附加源項，下標v、E、m分別代表動量、能量、質量；F為除SG內構件引入的附加力外的所有混合物所受體積力，如重力；keff為有效熱傳導系數；vdr,k為第k相的漂移速度。

動量、能量、空泡份額方程中的源項定義如下。動量方程中的源項表示由于控制體內SG構件的存在導致的附加阻力，包括下降段阻力、管束阻力、汽水分離器阻力、支承板阻力，其中，前3個阻力均為分布阻力，以分布阻力的形式添加到控制體網格，支承板阻力為局部阻力，以多孔階躍模型的方式添加到支承板所在位置的網格界面上。對順流和橫掠阻力分別采用管束外順流[14]和橫掠[15]阻力經驗關系式計算。支承板、下降段和汽水分離器的阻力，根據SG設計參數給定[13]。

能量方程的源項，即一次側向二次側釋熱，被簡化為一維分布計算。將管束等效為1根傳熱管，二次側流場也沿管束劃分一維控制體，對于直管段冷、熱側分別劃分11個控制體，彎管區域冷、熱側各劃分1個控制體，共計24個控制體。U型管內為單相液體對流，采用Dittus-Boelter公式[16]計算；管外換熱包括單相對流換熱、過冷沸騰換熱及飽和沸騰換熱，對于單相對流換熱，由于目前的經驗關系式多為管外順流和橫掠關系式，無傾斜沖刷管束的換熱關系式，因此將流速分解為順流和橫掠流速，分別計算順流換熱系數[17-18]和橫掠換熱系數[19]，將這兩者的加權和[13]作為實際的斜掠換熱系數。在計算中，給定一次側入口的質量流速和溫度，即可計算一、二次側間的換熱。

空泡份額方程中的源項表示氣相質量源項，其值等于蒸發率與冷凝率之差，對于穩態計算忽略SG內的非穩態冷凝現象，氣相質量源項等于液相蒸發率。

對于湍流模型，已公開發表的模型多針對純流體區域的精細網格提出[20]，對于多孔介質內的湍流模型，研究者在多孔介質粗網格內對已有湍流模型(如k-ε)進行積分得到適用于多孔介質的湍流模型[9-11]，但這些模型僅用于單相流動，并無針對兩相流動的湍流模型。本文采用文獻[21-22]提出的針對管束間兩相流動的零方程湍流模型計算湍流黏性及有效黏性，該模型已在ATHOS程序中驗證[13,23]。

以上模型構成完備的控制方程組，將這些方程在ANSYS FLUENT求解器中求解可獲得三維流場解，本文所采用的求解方法為COUPLE算法。

動員各方面力量搞好村莊環境宣傳，把全體村民發動起來，讓他們都真正理解支持，是每個村民都能做到愛村、護村，并行動起來，投入到美麗村莊建設上來，形成合力勢頭和建設發展力。

2 網格劃分及邊界條件

本文所選取的計算區域從下降段入口向下到圍板缺口，然后折流進入管束區，最終進入汽水分離器，計算區域示意圖如圖1a所示。為簡化網格劃分，采用方形通道代替圓筒形的汽水分離器筒體，保持筒體流通面積和中心位置不變。所有網格均為六面體網格，在支承板位置添加內部邊界面作為多孔階躍邊界。網格劃分如圖1b所示，網格質量大于0.5。

圖1 計算區域(a)及網格劃分(b)示意圖

網格進口邊界設置為速度進口，出口邊界設置為壓力出口，考慮到多孔介質網格尺寸較大，壁面無滑移邊界不適用，因此本文采用滑移壁面邊界條件，將壁面摩擦阻力分布添加到壁面網格內。入口溫度采用疏水和新給水的焓計算，出口壓力設置為SG頂部腔室壓力。

3 結果分析

采用FLUENT計算二次側流場、同時耦合一次側一維流場，得到二次側三維流場分布、一次側一維流場分布及一、二次側間耦合換熱系數分布。FLUENT計算值與設計值對比列于表1。

表1 FLUENT計算值與設計值的對比

3.1 一、二次側溫度分析

圖2為SG二次側流體溫度分布，圖2a中右側(x>0)為熱側、左側為冷側。由圖2a可看出，熱側流體在進入管束區后很快達到飽和溫度，而冷側流體緩慢達到飽和，冷側流體流過SG內1/2高度時才完全達到飽和，達到飽和后流體溫度保持不變。由圖2b的橫切面溫度分布可看出，由于流動和幾何的對稱性，流體溫度關于y=0平面對稱。

圖3為豎直線和水平線上的溫度分布。由圖3可看出，隨著軸向位置的升高冷卻劑溫度逐漸升高，且熱側(30、31曲線)流體升高速度遠較冷側(18曲線)的大，冷側流體在高度上升到7 m時才達到飽和，而熱側流體在管板附近即達到飽和。由對稱截面上水平線上的溫度分布可看出，當位于較低位置(z=0.2,1.0,3.0 m)時流體溫度在徑向分布很不均勻，徑向中心位置流體溫度較高，靠近外圍區域流體溫度較低，隨著流動的進行徑向外圍區域流體溫度升高，且熱側溫度高于冷側溫度。

圖4為流體及傳熱管溫度的一維分布。由圖4可看出，一次側流體溫度隨流動的進行逐漸降低，且隨沿管束流程的增加流體溫度降低速度減緩。熱側的二次側流體溫度沿軸向高度的增大很快達到飽和溫度，之后維持在飽和溫度544.7 K，冷側流體溫度沿軸向高度的增大緩慢達到飽和溫度。傳熱管溫度沿管束方向先增大后減小、之后又緩慢增大。這是由于在二次側流體溫度低于飽和溫度時，平均后的管外換熱屬于單相對流換熱，換熱系數較小；當流體溫度達到飽和后，沸騰換熱系數遠大于單相對流換熱系數。一、二次側的換熱系數如圖5所示，由于物性的變化，一次側單相對流換熱系數逐漸降低。二次側換熱系數在進入沸騰換熱區域后迅速增大，在沸騰區域換熱系數緩慢降低，直到再次進入單相對流區域又迅速降低。由于一、二次側換熱系數的突變，導致傳熱管內、外壁面溫度在單相對流區域升高。

a——對稱截面；b——不同高度水平截面

a——沿汽水分離器中心豎直線；b——沿對稱截面上水平線

圖4 流體及傳熱管溫度的一維分布

3.2 二次側流動含氣率分析

SG二次側流體的流動含氣率分布如圖6所示。由圖6可看出，在管束進口區域冷、熱側含氣率差別較大，熱側流動含氣率小于0.01的區域很小；而冷側該區域占較大面積。隨著流體沿管束上升，流體不斷蒸發，氣泡在橫向擴散，冷、熱側含氣率逐漸均衡，且熱側靠近中心處含氣率逐漸變為最大。

圖5 換熱系數分布

圖7為豎直線和水平線上的流動含氣率分布。流動含氣率隨軸向位置的增大而增大，熱側含氣率遠大于冷側的。同一高度上含氣率的最大值出現在熱側靠中心區域，含氣率的峰值隨軸向高度的增大逐漸向中心線移動。

a——對稱截面；b——不同高度水平截面

a——沿汽水分離器中心豎直線；b——沿對稱截面上水平線

圖8為汽水分離器內流動含氣率的分布。由圖8可看出，對于熱側最內層汽水分離器，其含氣率很高，達0.62，而對于冷側外圍的汽水分離器，含氣率僅0.05，不同的汽水分離器內混合物流動情況相差很大。

圖8 汽水分離器內流動含氣率分布

3.3 二次側壓力分析

圖9 對稱截面壓力分布

圖9為對稱截面上的壓力分布，壓力選取出口截面為參考面，該截面壓力設為0 Pa，圖中標注壓力為相對于參考壓力的相對壓力。由圖9可看出，由重位壓降、下降腔室壁面摩擦壓降和加速壓降的共同作用，下降段的壓力沿下降段向下逐漸增大；由于流體轉向，流體在圍板缺口區域壓力陡增，之后進入管束區域，流體壓力隨高度方向逐漸減小。由圖還可看出，同一高度截面上熱側壓力總體稍大于冷側壓力。圖10為壓力沿豎直線上的分布，壓力沿高度方向出現明顯的階躍變化，這是由于支承板的存在，引入局部阻力，導致流體經過支承板時，壓力突降；支承板導致壓降隨軸線高度的增大而增大，這是由于隨高度的增加，兩相流體含氣率增大，兩相摩擦壓降倍增因子增大，從而在相同的阻力系數下，支承板位置越高，引入的壓降越大。

圖10 壓力沿汽水分離器中心線的分布

3.4 二次側速度分析

圖11為對稱截面速度矢量圖。由圖11可看出，熱側流體速度明顯大于冷側的，最大速度約8 m/s。除圍板缺口折流部分和彎管區域，流體基本沿管束方向流動，橫流較小。管束區域流體速度矢量圖如圖12所示。

圖11 對稱截面速度矢量圖

豎直線和水平線上的速度分布如圖13所示。由圖13可看出，對于熱側流體，由于蒸發的作用，在管束區流體速度沿高度方向迅速上升；對于冷側流體，高度小于5 m時流體速度有輕微下降趨勢，這是由于冷側單相流體在流動過程中逐漸均勻化，5 m后由于蒸發的作用，流體速度開始上升。當流體流出管束區時，由于流通面積的增大，流體流速迅速降低，之后流入汽水分離器，由于流通面積減小，流體速度增大。由于冷、熱側加熱不均勻，流速在徑向上的分布也極不均勻。對于缺口高度的流體(z=0.2 m)，由于流體的折流，在外圍(即圍板缺口區域)流體流速較高，中心流速較低；隨高度的增大，熱側流體開始蒸發，流速逐漸升高，流速峰值出現在中心偏熱側區域。

圖12 對稱截面彎管部分速度矢量圖

a——沿汽水分離器中心豎直線；b——沿對稱截面上水平線

4 結語

本文采用多孔介質模型，對壓水堆SG二次側三維兩相流場進行數值模擬，同時耦合一、二次側一維換熱，得到如下結論：

1) 采用多孔介質模型計算SG二次側流場，相較于采用常規的CFD網格方法，可大幅節省計算量，并得到合理的數值解；

2) SG內冷熱側流體由于受熱不均，導致兩側流體流動狀態差別較大；

3) 汽水分離器內流動含氣率分布極不均勻，其值介于0.05和0.62之間，此分布可為汽水分離器設計提供進口邊界條件；

4) 通過計算得到彎管區速度場分布，可為U型管流致振動分析提供輸入參數。

參考文獻：

[1] PATANKAR S, SPALDING D. Heat exchanger design theory source book[M]. Washington D. C.: Scripta Book Co., 1974.

[2] SHA W, YANG C, KAO T, et al. Multidimensional numerical modeling of heat exchangers[J]. J Heat Transfer, 1982, 104(3): 417-425.

[3] PRITHIVIRAJ M, ANDREWS M. Three dimensional numerical simulation of shell-and-tube heat exchangers, Part Ⅰ: Foundation and fluid mechanics[J]. Numerical Heat Transfer, Part A: Applications, 1998, 33(8): 799-816.

[4] PRITHIVIRAJ M, ANDREWS M. Three-dimensional numerical simulation of shell-and-tube heat exchangers, Part Ⅱ: Heat transfer[J]. Numerical Heat Transfer, Part A: Applications, 1998, 33(8): 817-828.

[5] FERNG Y. 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.

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

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

[8] GETACHEW D, MINKOWYCZ W, LAGE J. A modified form of theκ-εmodel for turbulent flows of an incompressible fluid in porous media[J]. Int J Heat Mass Transfer, 2000, 43(16): 2 909-2 915.

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

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

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

[12] MASUOKA T, TAKATSU Y. Turbulence model for flow through porous media[J]. Int J Heat Mass Transfer, 1996, 39(13): 2 803-2 809.

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

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

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

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

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

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

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

[20] 陶文銓. 數值傳熱學[M]. 2版. 西安：西安交通大學出版社,2001.

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

[22] SCHLICHTING H, GERSTEN K. Boundary-layer theory[M]. 8th ed. Germany: Springer, 2000.

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