叢騰龍，田文喜，秋穗正，蘇光輝，謝永誠(chéng)，姚彥貴

(1.西安交通大學(xué) 能源與動(dòng)力工程學(xué)院，陜西 西安 710049；2.上海核工程研究設(shè)計(jì)院，上海 200233)

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

最早將多孔介質(zhì)模型應(yīng)用于換熱器的是英國(guó)學(xué)者Patankar等[1]，他們于20世紀(jì)70年代就采用此模型對(duì)換熱器殼側(cè)進(jìn)行了數(shù)值模擬。之后，Sha等[2]又在文獻(xiàn)[1]的基礎(chǔ)上模擬了SG和反應(yīng)堆堆芯中冷卻劑的流動(dòng)，Prithiviraj等[3-4]模擬了三維換熱器中的流動(dòng)。Ferng等[5-7]采用CFX軟件多孔介質(zhì)模型對(duì)SG管束區(qū)的二次側(cè)流場(chǎng)進(jìn)行了分析。

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

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

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

質(zhì)量守恒方程：

β·(βρmvm)=0

(1)

β·(βρmvmvm)=-β

(2)

能量守恒方程：

β·

(3)

空泡份額方程：

β·(βαgρgvm)=

(4)

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

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

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

空泡份額方程中的源項(xiàng)表示氣相質(zhì)量源項(xiàng)，其值等于蒸發(fā)率與冷凝率之差，對(duì)于穩(wěn)態(tài)計(jì)算忽略SG內(nèi)的非穩(wěn)態(tài)冷凝現(xiàn)象，氣相質(zhì)量源項(xiàng)等于液相蒸發(fā)率。

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

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

動(dòng)員各方面力量搞好村莊環(huán)境宣傳，把全體村民發(fā)動(dòng)起來(lái)，讓他們都真正理解支持，是每個(gè)村民都能做到愛(ài)村、護(hù)村，并行動(dòng)起來(lái)，投入到美麗村莊建設(shè)上來(lái)，形成合力勢(shì)頭和建設(shè)發(fā)展力。

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

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

圖1 計(jì)算區(qū)域(a)及網(wǎng)格劃分(b)示意圖

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

3 結(jié)果分析

采用FLUENT計(jì)算二次側(cè)流場(chǎng)、同時(shí)耦合一次側(cè)一維流場(chǎng)，得到二次側(cè)三維流場(chǎng)分布、一次側(cè)一維流場(chǎng)分布及一、二次側(cè)間耦合換熱系數(shù)分布。FLUENT計(jì)算值與設(shè)計(jì)值對(duì)比列于表1。

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

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

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

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

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

a——對(duì)稱(chēng)截面；b——不同高度水平截面

a——沿汽水分離器中心豎直線；b——沿對(duì)稱(chēng)截面上水平線

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

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

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

圖5 換熱系數(shù)分布

圖7為豎直線和水平線上的流動(dòng)含氣率分布。流動(dòng)含氣率隨軸向位置的增大而增大，熱側(cè)含氣率遠(yuǎn)大于冷側(cè)的。同一高度上含氣率的最大值出現(xiàn)在熱側(cè)靠中心區(qū)域，含氣率的峰值隨軸向高度的增大逐漸向中心線移動(dòng)。

a——對(duì)稱(chēng)截面；b——不同高度水平截面

a——沿汽水分離器中心豎直線；b——沿對(duì)稱(chēng)截面上水平線

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

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

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

圖9 對(duì)稱(chēng)截面壓力分布

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

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

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

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

圖11 對(duì)稱(chēng)截面速度矢量圖

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

圖12 對(duì)稱(chēng)截面彎管部分速度矢量圖

a——沿汽水分離器中心豎直線；b——沿對(duì)稱(chēng)截面上水平線

4 結(jié)語(yǔ)

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

1) 采用多孔介質(zhì)模型計(jì)算SG二次側(cè)流場(chǎng)，相較于采用常規(guī)的CFD網(wǎng)格方法，可大幅節(jié)省計(jì)算量，并得到合理的數(shù)值解；

2) SG內(nèi)冷熱側(cè)流體由于受熱不均，導(dǎo)致兩側(cè)流體流動(dòng)狀態(tài)差別較大；

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

4) 通過(guò)計(jì)算得到彎管區(qū)速度場(chǎng)分布，可為U型管流致振動(dòng)分析提供輸入?yún)?shù)。

參考文獻(xiàn)：

[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] 陶文銓. 數(shù)值傳熱學(xué)[M]. 2版. 西安：西安交通大學(xué)出版社,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.