基于流固共軛傳熱的兩相流動數(shù)值模擬及燃料子通道CHF預(yù)測研究

馮琳娜，李 權(quán),*，黃永忠，劉振海，齊飛鵬，M. Avramova

(1.中國核動力研究設(shè)計院 核反應(yīng)堆系統(tǒng)設(shè)計技術(shù)重點(diǎn)實(shí)驗室，四川 成都 610213；2.北卡羅來納州立大學(xué)，美國 羅利 27695)

對于壓水堆燃料組件設(shè)計，臨界熱流密度(CHF)是評價其熱工水力性能的重要參數(shù)之一，一旦燃料棒的熱流密度超過CHF，就有可能發(fā)生燃料棒燒毀，進(jìn)而導(dǎo)致放射性物質(zhì)釋放的風(fēng)險。因此，對于燃料組件設(shè)計，需掌握其CHF值[1]。目前CHF值主要通過試驗的方法獲取，但試驗成本高、周期長，研究者一直在尋找替代試驗的數(shù)值預(yù)測方法。

對CHF預(yù)測的研究由來已久，如經(jīng)驗關(guān)系式、理論模型(汽泡雍塞模型等)、CHF查詢表等，但經(jīng)驗關(guān)系式由試驗數(shù)據(jù)擬合得到，只適用于特定的燃料組件，對于新的燃料組件設(shè)計則無法適用；理論模型、CHF查詢表等用于燃料組件復(fù)雜結(jié)構(gòu)則需引入一些修正因子，無法準(zhǔn)確評估燃料組件細(xì)節(jié)特征的影響[2]。因此隨著近年來CFD技術(shù)的興起，研究者一直在嘗試采用兩相CFD方法進(jìn)行CHF預(yù)測。Vyskocil等[3]利用NEPTUNE_CFD程序，采用兩相CFD分析的方法，認(rèn)為空泡份額超過0.8時CHF發(fā)生，并對豎直加熱的圓管CHF值進(jìn)行了預(yù)測。李權(quán)等[4-5]采用沸騰曲線的方法，對均勻加熱圓管和非均勻加熱圓管的CHF值進(jìn)行了預(yù)測，并探究了不同工況的影響。李松蔚等[6]采用不考慮相間質(zhì)量傳遞的兩相CFD分析的方法，分析了帶7個格架的5×5棒束通道兩相流動特性，并通過格架下游的空泡份額分布定性地判斷CHF特性。但上述數(shù)值模擬都只針對流體區(qū)域進(jìn)行求解，熱工邊界較為簡單，未能考慮流固共軛傳熱引起的軸向非均勻傳熱特性的影響，也無法準(zhǔn)確捕捉CHF發(fā)生時固體表面的溫度變化。而陳廣亮等[7]的研究表明，采用流固共軛傳熱可更加準(zhǔn)確地預(yù)測周向的非均勻傳熱特性。

因此，本文基于流固共軛傳熱和兩相CFD分析的方法，建立相應(yīng)的數(shù)值模型，實(shí)現(xiàn)流固共軛傳熱和過冷沸騰兩相流動的耦合求解，并在此基礎(chǔ)上建立燃料子通道CHF的數(shù)值預(yù)測方法。

1 數(shù)值模型

1.1 固體導(dǎo)熱模型

采用穩(wěn)態(tài)計算，固體導(dǎo)熱方程如下：

(1)

(2)

式中：λs為固體的熱導(dǎo)率；Ts為固體溫度；Qv為體熱源。

1.2 流固共軛傳熱數(shù)值模型

流固共軛傳熱主要通過流體和固體的交界面進(jìn)行傳遞，交界面上固體和流體的溫度、熱流密度保持連續(xù)。

Ts|surface=Tl|surface

(3)

qs|surface=ql|surface

(4)

式中：Ts|surface和Tl|surface分別為交界面上固體和流體的溫度；qs|surface和ql|surface分別為交界面上固體域和流體域的熱流密度。

1.3 過冷沸騰兩相流動數(shù)值模擬

對過冷沸騰兩相流動的CFD模擬，以歐拉兩流體模型為框架，結(jié)合壁面沸騰模型進(jìn)行。

1) 歐拉兩流體模型

歐拉兩流體模型是將每一相均看成充滿整個流動空間的連續(xù)介質(zhì)，用空泡份額來表征每相所占據(jù)的體積份額，對兩相分別建立質(zhì)量、動量和能量守恒方程，通過兩相間的界面?zhèn)鬟f模型使方程組封閉。建立的守恒方程為：

(5)

(6)

(7)

2) 兩相間的動量傳遞

兩相間的動量傳遞主要通過界面力實(shí)現(xiàn)，本文主要考慮曳力和湍流耗散力。

Mk=FD,k+FTD,k

(8)

式中，F(xiàn)D,k、FTD,k分別為由曳力、湍流耗散力引起的兩相間的動量傳遞。

曳力FD通常包含兩部分：時均值和脈動值。時均部分的曳力計算關(guān)系式為：

(9)

式中：ur=ul-ug為兩相間的相對速度；CD為曳力系數(shù)，一般由Tomiyama關(guān)系式[8]計算；指數(shù)n取4，用于考慮高空泡份額的效應(yīng)，也稱汽泡云集效應(yīng)[9]。

曳力的脈動部分，也稱為湍流耗散力，表征由液相湍流引起的動量傳遞。其值由下式計算：

(10)

3) 汽泡直徑分布

汽泡脫離壁面進(jìn)入主流區(qū)后，會與主流區(qū)的液體進(jìn)行熱交換，汽泡直徑會發(fā)生變化。由于汽泡在主流區(qū)主要發(fā)生冷凝，因此通常汽泡直徑為主流液體過冷度的線性關(guān)系式：

(11)

式中：dB為汽泡在主流區(qū)的直徑分布；Tsub為液體的過冷度。本文參考文獻(xiàn)[10]，取Tsub,1=13.5 K，dB1=0.1 mm；Tsub,2=5 K，dB2=2 mm。

4) 壁面沸騰模型

對過冷沸騰的模擬，還需加熱壁面的源項。為考慮高空泡份額及臨界熱流密度的影響，本文采用拓展后的壁面沸騰模型，將壁面熱流密度分成4部分，機(jī)理示意圖如圖1所示。

圖1 熱流密度分配示意圖[7]Fig.1 Diagram of heat flux distribution[7]

(12)

式中，qw=qs|surface為加熱壁面的熱流密度；qc、qq、qe、qg分別為液體單相對流傳熱、激冷對流傳熱、汽化潛熱及氣體單相對流傳熱部分的熱流密度，其詳細(xì)表達(dá)式參見文獻(xiàn)[9]；Kdry為壁面燒干面積比，由式(13)計算。

(13)

(14)

式中：αdry為臨界空泡份額；αδ為近壁面空泡份額。本文參考Weisman & Pei的汽泡雍塞理論[11]，αdry取值為0.82，αδ則取距壁面最近的網(wǎng)格柵元的值。

式(13)中各項的表達(dá)式參見文獻(xiàn)[5]，與文獻(xiàn)[5]不同的是，本文活化核心密度由Lemmert-Chawla模型[12]計算：

(15)

式中：C=210；m=1.805；Ts,sup=Ts-Tsat為加熱壁面的過熱度，Tsat為液體飽和溫度。

汽泡脫離直徑由Tolubinsky模型[13]計算：

Dd=d0exp(-Tsub/45)

(16)

式中：Dd為汽泡脫離直徑；d0為參考直徑。本文主要模擬高壓工況，高壓下汽泡脫離直徑更小，因此參考Krepper等[10]的研究取d0=0.6 mm。

上述數(shù)值模型均應(yīng)用于軟件STAR-CCM+9.04[14]。

2 燃料棒典型柵元共軛傳熱兩相流動分析

2.1 數(shù)值模型建立

沸水堆中冷卻劑處于沸騰兩相流動狀態(tài)，可為本文建立的數(shù)值模型提供驗證。選擇典型的BWR6沸水堆燃料棒[14]，該燃料棒長度為3.7 m，燃料芯塊直徑為10.68 mm，燃料芯塊的功率密度為226.4 MW/m3，包殼壁厚為0.81 mm，冷卻劑運(yùn)行壓力為7.17 MPa，本文取該燃料棒前1 m進(jìn)行模擬，該階段處于過冷沸騰兩相流動狀態(tài)。詳細(xì)的幾何模型如圖2所示。

圖2 幾何模型及網(wǎng)格示意圖Fig.2 Diagram of geometrical model and mesh

在數(shù)值模擬中，采用二維旋轉(zhuǎn)對稱模型，對該幾何模型進(jìn)行建模。在燃料芯塊與包殼之間，通過增加熱阻的方式來模擬燃料棒中預(yù)填充氦氣的影響。

采用結(jié)構(gòu)化網(wǎng)格進(jìn)行網(wǎng)格劃分，經(jīng)過網(wǎng)格敏感性分析，確定總的網(wǎng)格數(shù)量為13 600，局部網(wǎng)格示意圖如圖2所示。

2.2 模型驗證

基于上述網(wǎng)格和所建立的數(shù)值模型，進(jìn)行燃料芯塊、包殼、冷卻劑的共軛傳熱兩相模擬，計算得到靠近出口處的芯塊溫度、包殼溫度和冷卻劑溫度沿徑向的分布與理論值[14]的對比如圖3所示。從圖3可見，數(shù)值模擬得到的芯塊溫度、包殼溫度和流體溫度均與理論值符合良好，驗證了本文基于流固共軛傳熱的過冷沸騰兩相流動數(shù)值模型和分析方法的正確性。

圖3 芯塊、包殼及冷卻劑溫度沿徑向的變化Fig.3 Radial temperature distribution of fuel pellet, cladding and coolant

3 典型燃料子通道內(nèi)兩相流動及CHF數(shù)值預(yù)測

3.1 幾何建模及網(wǎng)格劃分

在典型柵元的基礎(chǔ)上，探索燃料子通道內(nèi)兩相流動和CHF的預(yù)測方法，以國際性基準(zhǔn)驗證實(shí)驗PSBT[15]中子通道S1的試驗數(shù)據(jù)進(jìn)行模擬。子通道S1包含4塊1/4加熱片，其柵元距離為12.6 mm，結(jié)構(gòu)示意圖如圖4a所示，加熱元件采用鉻鎳鐵合金，直徑為9.5 mm，在子通道外采用礬土進(jìn)行絕緣。加熱段長度均為1.555 m，為均勻加熱。試驗測量了加熱段入口1.4 m處截面的平均空泡份額。

圖4 子通道S1及網(wǎng)格劃分尺寸示意圖Fig.4 Schematic diagram of subchannel S1 and partial grids

由于幾何模型的對稱性，采用1/8子通道模型進(jìn)行建模并施加對稱邊界條件，同時對加熱元件和流體域進(jìn)行建模。軸向上分別增加20 mm的入口段和出口段，以排除入口和出口效應(yīng)的影響。采用結(jié)構(gòu)化網(wǎng)格進(jìn)行數(shù)值模擬，通過網(wǎng)格敏感性分析，最終選擇的網(wǎng)格示意圖如圖4b所示。

3.2 兩相流動數(shù)值模擬

基于上述網(wǎng)格和數(shù)值模型，對表1中14個工況進(jìn)行數(shù)值模擬，計算得到距加熱段1.4 m高度處的截面平均空泡份額計算值與試驗值的對比如圖5所示。從圖5可見，計算值與試驗值總體符合良好，表明本文模型可用于燃料子通道過冷沸騰兩相流動模擬。

表1 子通道兩相流動模擬工況Table 1 Test condition of boiling two-phase flow in subchannel

3.3 典型子通道CHF預(yù)測

在兩相流動模擬的基礎(chǔ)上，對CHF的預(yù)測進(jìn)行研究。以1.2221工況為例，采用準(zhǔn)瞬態(tài)的方法，保持入口條件不變，逐步增加加熱功率，并監(jiān)測加熱元件的最大壁面溫度，計算得到的最大壁面溫度隨熱流密度(加熱功率轉(zhuǎn)換得到)的變化如圖6所示。從圖6可見，在進(jìn)入到過冷沸騰區(qū)后，剛開始時，壁面最大溫度隨熱流密度的提升變化幅度較小。但隨著加熱功率的提高，當(dāng)達(dá)到一定的熱流密度后，壁面溫度升高的趨勢明顯增大，有溫度飛升的趨勢，此現(xiàn)象與臨界熱流密度試驗中出現(xiàn)壁面溫度飛升的現(xiàn)象非常近似，因此，本文將變化曲線中壁面溫度出現(xiàn)飛升作為臨界熱流密度發(fā)生的判據(jù)。則對于1.2221工況，CHF預(yù)測值為2 400 kW/m2。

圖5 兩相流動模擬工況計算值與試驗值的對比Fig.5 Comparison between calculated and test values of two-phase flow simulation

圖6 1.2221工況的最大壁面溫度隨熱流密度的變化Fig.6 Change of maximum wall temperature with heat flux in 1.2221 condition

在PSBT子通道實(shí)驗中，并無測量相關(guān)的CHF實(shí)驗數(shù)據(jù)。本文借助2006年的CHF查詢表[16]，計算得到CHF的實(shí)驗值。對于1.2221工況，CHF查詢表值為2 088 kW/m2，預(yù)測誤差為14.9%。造成CHF預(yù)測值較查詢表值大的因素有兩個：1) 熱流密度為2 400 kW/m2時，出口的局部含汽率達(dá)0.075，根據(jù)文獻(xiàn)[5]研究結(jié)果，在較高含汽率下，預(yù)測值容易較實(shí)驗值高；2) 熱流密度2 400 kW/m2時，壁面最大空泡份額已超過0.82，此時計算得到的局部含汽率較實(shí)際CHF值偏高，導(dǎo)致通過CHF查詢表得到的CHF值偏低。

基于壁面溫度飛升的方法，對表1中1.2221、1.1221、1.2211、1.2231、1.2233、1.2234、1.2423、1.2121、1.3221共9個工況的CHF值進(jìn)行了預(yù)測，CHF預(yù)測值與查詢表值的對比如圖7所示。

圖7 子通道CHF預(yù)測值與查詢表值的對比Fig.7 Comparison of predicted CHF with look-up table data

由于上述兩個因素，對于大部分工況，預(yù)測值均較CHF查詢表值高，9個工況的平均預(yù)測誤差為15%。但對于1.2423和1.3221工況，預(yù)測值較CHF查詢表值低。對于1.2423工況，流量為1 369.4 kg/(m2·s)，明顯小于其他工況；而1.3221工況的壓力為12.3 MPa，也小于其他工況；文獻(xiàn)[5]發(fā)現(xiàn)，采用同一套數(shù)值模型對CHF進(jìn)行預(yù)測，隨著壓力和流量的增加，預(yù)測值與實(shí)驗值的比值逐漸增加，即在相對低壓、低流量下，預(yù)測值較實(shí)驗值低，因此導(dǎo)致這兩個工況的預(yù)測值較實(shí)驗值低。本文模擬的趨勢與文獻(xiàn)[5]相符合。

以上分析表明，基于流固共軛傳熱兩相流動模擬構(gòu)建的CHF預(yù)測方法，可適用于子通道內(nèi)CHF的數(shù)值預(yù)測，但本文CHF實(shí)驗值的獲取方式仍需進(jìn)一步提高，以更加準(zhǔn)確地評估本文預(yù)測方法的精度。

4 結(jié)論

本文基于流固共軛傳熱和兩相CFD分析的方法，建立了考慮共軛傳熱的過冷沸騰兩相流動數(shù)值模型，并通過典型沸水堆燃料棒的分析，驗證了數(shù)值模型的正確性。在此基礎(chǔ)上，對燃料子通道的兩相流動進(jìn)行了模擬，并采用準(zhǔn)穩(wěn)態(tài)求解的方法，基于加熱壁面溫度飛升，建立了CHF預(yù)測方法。采用該方法，對不同工況燃料棒束子通道內(nèi)的CHF進(jìn)行數(shù)值模擬，驗證了CHF預(yù)測方法的正確性。本文基于共軛傳熱兩相流動的CHF預(yù)測方法可為后續(xù)燃料組件和格架等部件的設(shè)計優(yōu)化提供參考。