吳宗蕓，劉天才，吳明宇

(中國原子能科學(xué)研究院，北京 102413)

由于子通道分析方法能詳細(xì)地考慮反應(yīng)堆組件內(nèi)每個通道的質(zhì)量、能量、動量守恒過程并且相對于CFD方法具有較高的計(jì)算效率，因此子通道分析方法被廣泛用于鈉冷快堆的設(shè)計(jì)。然而，由于鈉冷快堆普遍使用繞絲對燃料棒進(jìn)行定位，對鈉冷快堆組件進(jìn)行子通道分析時，需對繞絲帶來的交混效應(yīng)進(jìn)行準(zhǔn)確模擬以更準(zhǔn)確計(jì)算冷卻劑的溫度場分布。在子通道分析程序中常使用的繞絲模型可分為3類：強(qiáng)迫橫流模型、分布式阻力模型和雙區(qū)域模型。在COBRA[1]系列子通道程序中，繞絲對組件內(nèi)溫度場以及流場帶來的影響使用強(qiáng)迫橫流模型來處理，這種模型認(rèn)為繞絲在繞過燃料棒間的間隙時，會帶來沿著繞絲繞向的橫向流動。在MATRA-LMR[2]、SACOS-PB[3]、ATHAS-LMR[4]等子通道分析程序中，使用分布式阻力模型[4]來模擬繞絲帶來的效應(yīng)，這類模型通過考慮繞絲對軸向流動以及間隙橫流的流動阻力來模擬繞絲對流場的影響。相對于強(qiáng)迫橫流模型，分布式阻力模型更適用于低流量工況。但這兩種模型在組件徑向功率傾斜分布時，計(jì)算得到的結(jié)果與實(shí)驗(yàn)值相差較大[5]。第3類在子通道分析中使用的繞絲模型是雙區(qū)域模型[6]，這種模型被用在專門針對鈉冷快堆開發(fā)的子通道分析程序ENERGY和SLTHEN中，該模型可較為詳細(xì)地考慮子通道間的各種交混機(jī)理并且相對于前兩種模型能較為準(zhǔn)確預(yù)測組件非均勻功率分布下的冷卻劑溫度場[7]。但由于ENERGY和SLTHEN程序中使用簡化后的子通道守恒方程，ENERGY和SLTHEN程序使用理論計(jì)算來得到割流系數(shù)以及內(nèi)部通道和壁面通道的速度，并且程序內(nèi)部的模型只有內(nèi)部區(qū)域和外部區(qū)域兩個速度場，每個區(qū)域有相同的流速，而在實(shí)際情況中每個通道內(nèi)的流速往往是不同的。因此ENERGY和SLTHEN程序計(jì)算得到的結(jié)果與一些重要的實(shí)驗(yàn)數(shù)據(jù)有較大偏差[7]。

除上述程序外，針對液態(tài)金屬冷卻快堆子通道分析，近年來還開發(fā)了一系列其他的子通道分析程序，如KMC-Sub[8]、SUBAC[9]和COBRA-LM[10]，但這幾款程序使用集總的交混系數(shù)來考慮繞絲的湍流交混，模型較為粗糙；ANTEO子通道分析程序[11]中也使用了雙區(qū)域模型來考慮繞絲的交混效應(yīng)，其對WARD實(shí)驗(yàn)的驗(yàn)證結(jié)果在高流量工況下明顯好于二次開發(fā)后的COBRA，但該程序在低雷諾數(shù)工況下與實(shí)驗(yàn)值有較大的偏差[11]；一些其他開發(fā)的子通道分析程序的文獻(xiàn)[12-14]中未深入介紹使用的繞絲交混模型。在這些程序公開發(fā)表的驗(yàn)證工作中只有極少數(shù)對功率偏斜工況進(jìn)行了驗(yàn)證，而且國內(nèi)關(guān)于雙區(qū)域模型的子通道分析程序的開發(fā)研究在公開文獻(xiàn)中還較為少見。

本文使用雙區(qū)域模型來模擬繞絲帶來的效應(yīng)，建立相應(yīng)的子通道守恒方程來考慮每個通道內(nèi)的質(zhì)量、能量、動量守恒，并結(jié)合液態(tài)金屬鈉的物性、對流換熱與流動阻力模型開發(fā)專門針對鈉冷快堆組件熱工水力分析的子通道程序SPLICA(sub-channel program for liquid metal cooled assembly)。

1 帶繞絲棒束組件的雙區(qū)域模型子通道模型

在雙區(qū)域模型中，整個組件被分為兩個區(qū)域：由內(nèi)部通道組成的內(nèi)部區(qū)域以及由邊通道和角通道組成的外部區(qū)域。在該模型中，繞絲在內(nèi)部區(qū)域和外部區(qū)域的交混效應(yīng)采用不同的方式考慮。在由內(nèi)部通道組成的內(nèi)部區(qū)域中，間隙兩端圍繞燃料棒螺旋的繞絲交替穿過間隙，導(dǎo)致了間隙上的間隙橫流強(qiáng)度的一次諧波近似為正弦函數(shù)[15]。雙區(qū)域模型示意圖如圖1所示，繞絲在內(nèi)部區(qū)域周期性的間隙橫流增強(qiáng)了通道間的擴(kuò)散作用，使子通道間的溫差與冷卻劑流量差減小。而在由邊通道和交通道的外部區(qū)域中，由于繞絲每次穿過燃料棒與組件盒之間的外圍間隙時均沿著同一方向，因此帶來了外圍間隙的單向?qū)α鹘换?。對于處于堆芯外圍的組件，燃料棒的功率在組件內(nèi)呈現(xiàn)明顯的傾斜非均勻分布，這時外圍間隙的單向流動對組件內(nèi)冷卻劑的溫度分布具有非常明顯的影響。利用這樣的基本物理模型，建立子通道分析雙區(qū)域模型的基本守恒方程與本構(gòu)模型。

圖1 雙區(qū)域模型示意圖Fig.1 Schematic diagram of two-region model

1.1 子通道守恒方程

子通道控制體質(zhì)量守恒方程：

(1)

能量守恒方程：

(2)

軸向動量守恒方程：

(3)

與COBRA這類子通道分析程序不同，本文開發(fā)的子通道分析程序不包含橫向動量守恒方程。由于在鈉冷快堆中，繞絲對間隙的橫流強(qiáng)度與方向起著決定性的影響，因此使用橫向動量守恒方程來求解得到的間隙橫流往往與實(shí)驗(yàn)值具有較大的偏差。在本文開發(fā)的子通道分析程序中，軸向動量守恒方程采用軸向等壓近似求解[15]。

1.2 繞絲棒束組件的流動壓降模型

為使子通道分析程序能適用于低流量工況下的計(jì)算，流動壓降模型需考慮層流區(qū)、層流-湍流過渡區(qū)、湍流區(qū)下的流動阻力。本文選擇詳細(xì)CT(Cheng-Todreas)模型[16]來計(jì)算繞絲棒束的流動壓降，該模型覆蓋了廣泛的雷諾數(shù)范圍，具有良好的適用性。該模型不同類型的通道流動阻力使用不同的公式表示。

內(nèi)部通道的阻力系數(shù)f1：

(4)

邊通道的阻力系數(shù)f2：

(5)

角通道的阻力系數(shù)f3：

(6)

式中：m為常數(shù)，層流流動工況下m=1.0，湍流流動工況下m=0.18；θ為繞絲與軸向方向的夾角；下標(biāo)i=1,2,3分別表示內(nèi)部通道、邊通道以及角通道。

當(dāng)棒束通道的平均雷諾數(shù)RebRebT時認(rèn)為流動處在湍流區(qū)。當(dāng)RebL

(7)

對于層流和湍流過渡區(qū)，流動阻力系數(shù)f使用插值計(jì)算：

f=fRebLψ1/3+fRebT(1-ψ)1/3

(8)

1.3 單相對流換熱模型

在本程序中單相對流換熱系數(shù)采用Mikityuk關(guān)系式[17]。Mikityuk在2009年對以前的對流換熱關(guān)系式進(jìn)行了對比分析，擬合了如下的關(guān)系式：

Nu=0.047(1-exp(-3.8(P/d-1)))·

(Pe0.77+250)

(9)

式中：S為燃料棒的間距；d為燃料棒的直徑；Pe為佩特萊數(shù)。適用范圍：1.1≤P/d≤1.95，30≤Pe≤5 000。

1.4 繞絲棒束通道間的交混模型

本文選用Zhukov交混模型和Cheng-Todreas交混模型來計(jì)算通道間的能量交混系數(shù)和動量交混系數(shù)。交混模型是封閉雙區(qū)域子通道守恒方程組的關(guān)鍵，并對計(jì)算得到的結(jié)果的準(zhǔn)確性有至關(guān)重要的影響。

1)Zhukov繞絲棒束交混模型

Zhukov交混模型使用下面的公式計(jì)算交混系數(shù)。

(10)

式(10)的適用范圍為40≤Pe≤1 500,1.15≤S/d≤1.32,0.005≤Pr≤0.03，其中，Re為間隙相連兩個通道的平均雷諾數(shù)，Pr為平均普朗特?cái)?shù)。

Ψ(Re)=1-0.694exp(0.132×10-3Re)

(11)

(12)

式中，Δ為組件盒內(nèi)壁面與最外側(cè)燃料棒包殼之間的間隙寬度。

(13)

(1-exp(-40(S/d-1)))

(14)

2)Cheng-Todreas繞絲棒束交混模型

與Zhukov模型不同，Cheng-Todreas模型使用一個系數(shù)來考慮內(nèi)部通道的間隙湍流與繞絲交混引起的總體擴(kuò)散效應(yīng)，而非將分子-湍流交混和繞絲引起的交混單獨(dú)考慮。并且Cheng-Todreas模型依賴于棒束的流動是否處于湍流區(qū)，對于湍流和層流流動工況使用不同的公式計(jì)算，在層流和湍流過渡區(qū)時，Cheng-Todreas模型使用插值來計(jì)算交混系數(shù)。Cheng-Todreas模型[23]采用如下的公式來計(jì)算通道之間的交混：

(15)

式中：Gin為內(nèi)部通道的質(zhì)量流量；θ為繞絲與軸向方向的夾角；Bij為間隙的寬度。

對于組件的外圍間隙，其與1根燃料棒以及燃料組件盒壁相接，因此湍流帶來的擴(kuò)散效應(yīng)以及繞絲帶來的單向?qū)α鹘换煨鑶为?dú)考慮。

外圍間隙的湍流擴(kuò)散交混使用下式計(jì)算：

(16)

式中，Gside為邊通道的質(zhì)量流量。

外圍間隙的繞絲帶來的單向?qū)α鹘换焓褂孟率接?jì)算：

C1L=Cm1(Ar2/A′1)0.5tanθ

(17)

式中，系數(shù)Cm1、Cm2依賴于棒束的流動工況以及組件棒束的幾何。

對于湍流區(qū)(Reb>RebT)以及組件燃料棒數(shù)Nr≥19：

Cm1T=0.14(B/D)-0.5

Cm2T=0.75(H/D)0.3

(18)

對于湍流區(qū)以及組件燃料棒數(shù)Nr≤7：

Cm1T=0.1(B/D)-0.5

Cm2T=0.6(H/D)0.3

(19)

對于層流區(qū)(Reb

Cm1L=0.5(B/D)0.6Cm1T

Cm2L=0.5(B/D)0.6Cm2T

(20)

對于層流湍流過渡區(qū)(RebL

Cm1=Cm1L+(Cm1T-Cm1L)ψ2/3

Cm2=Cm2L+(Cm2T-Cm2L)ψ2/3

(21)

對于Cheng-Todreas交混模型，動量交混系數(shù)ηk,j與能量交混系數(shù)wk,j的取值一致。

2 基于雙區(qū)域模型的子通道分析程序的開發(fā)

本文利用面向?qū)ο蟮腃++語言完成子通道分析程序SPLICA的開發(fā)，程序中能量守恒方程在組件的每個軸向平面上使用SOR迭代法來隱式求解，動量守恒方程與質(zhì)量守恒方程聯(lián)立并采用軸向等壓近似方法求解。SPLICA程序可接受系統(tǒng)壓力、入口流量或組件棒束壓降、入口冷卻劑溫度或冷卻劑比焓作為其邊界條件。程序具有穩(wěn)態(tài)以及緩慢變化瞬態(tài)的計(jì)算功能，并且能計(jì)算燃料元件內(nèi)部的熱傳導(dǎo)。圖2為本文所開發(fā)的子通道分析程序的計(jì)算流程圖。程序使用逐層求解的方式，在每一層計(jì)算完成后再沿著軸向計(jì)算下一層，最后得到整個組件上的熱工水力工況參數(shù)。在每次外迭代中，程序求解燃料元件的導(dǎo)熱微分方程以更新燃料元件表面的熱流密度。并從組件的入口層沿著軸向計(jì)算到出口層，對每層冷卻劑控制體使用能量守恒方程求解得到冷卻劑的溫度，然后根據(jù)當(dāng)前層冷卻劑控制體的溫度與壓力更新當(dāng)前層冷卻劑控制體的物性，之后根據(jù)動量方程與質(zhì)量守恒方程，求解得到當(dāng)前層控制體的出口冷卻劑質(zhì)量流量。當(dāng)整個外迭代收斂時，輸出計(jì)算結(jié)果。

圖2 SPLICA程序計(jì)算流程圖Fig.2 Flow chart of SPLICA subchannel code

3 算例驗(yàn)證與測試

3.1 ORNL-FFM 19棒束實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證

FFM(fuel failure mockup)組件實(shí)驗(yàn)[24]是由美國橡樹林國家實(shí)驗(yàn)室(ORNL)進(jìn)行的。該實(shí)驗(yàn)裝置專門用于研究鈉冷快堆的堆芯熱工水力問題。FFM是一個高溫鈉冷卻實(shí)驗(yàn)裝置，測試段使用19根棒束模擬液態(tài)金屬冷卻反應(yīng)堆的堆芯組件，采用電加熱棒將熱量傳遞給冷卻劑鈉。該實(shí)驗(yàn)裝置的測試段子通道劃分與編號方案如圖3所示。FFM實(shí)驗(yàn)中液態(tài)金屬鈉的溫度最高可達(dá)650 ℃，棒束的最大線功率為33.3 kW/m，冷卻劑最高流量為12.2 kg/s、最低為0.05 kg/s。整個實(shí)驗(yàn)段的長度為1.016 m，其中非加熱入口段長305 mm，加入實(shí)驗(yàn)段長度為533.75 mm，非加熱出口段的長度為76.25 mm。實(shí)驗(yàn)段組件的參數(shù)列于表1。表2列出了101～105組高流量實(shí)驗(yàn)的入口溫度、流量以及加熱棒功率條件，整個測試段運(yùn)行在大氣壓下。

圖3 FFM 2A實(shí)驗(yàn)組件編號方案Fig.3 Sub-channel numbering scheme for FFM 2A assembly

表1 FFM 2A 19棒束實(shí)驗(yàn)測試段組件參數(shù)Table 1 FFM 2A 19-rod bundle experiment test section component parameter

表2 FFM 2A 高流量實(shí)驗(yàn)工況參數(shù)Table 2 Parameter of FFM 2A high flow rate experimental condition

圖4示出了對于FFM 2A高流量實(shí)驗(yàn)102～105組計(jì)算得到的無量綱溫度的計(jì)算值與實(shí)驗(yàn)值，子通道i的無量綱溫度Tnorm,i的定義為：

(22)

圖4 FFM 2A實(shí)驗(yàn)段出口無量綱溫度計(jì)算值與實(shí)驗(yàn)值對比Fig.4 Comparison between calculated and experimental values of dimensionless temperature at exit of FFM 2A experiment section

3.2 WARD 61棒束實(shí)驗(yàn)數(shù)據(jù)數(shù)據(jù)驗(yàn)證

WARD(westinghouse advanced reactors division)實(shí)驗(yàn)段的測試組件擁有61根加熱棒[25]。WARD實(shí)驗(yàn)測試段棒束總長度為265 cm，其中加熱段的長度為114.3 cm。冷卻劑鈉從底端入口流入，經(jīng)過加熱段受熱后從頂端流出。在加熱段的底端有一個長24.1 cm的非加熱入口段。加熱段的軸向功率分布為截?cái)嗟挠嘞曳植?，加熱段軸向線功率最大值與平均值的比值為1.40。實(shí)驗(yàn)測量了徑向上均勻的功率分布以及傾斜的功率分布下，加熱段以及加熱段下游的非加熱區(qū)共5個截面上的通道溫度。由于WARD實(shí)驗(yàn)測量了大量的層流工況以及層流-湍流過渡工況下的實(shí)驗(yàn)數(shù)據(jù)，因此WARD實(shí)驗(yàn)可用來校驗(yàn)子通道程序在低流量工況下的溫度分布計(jì)算結(jié)果。圖5示出了WARD 61棒束實(shí)驗(yàn)的子通道劃分與編號。在計(jì)算時，選取測量距離加熱段底端57.2、115.6、179.1 cm的測量數(shù)據(jù)和計(jì)算數(shù)據(jù)進(jìn)行比較。其中距加熱段底端57.2 cm的測量截面位于加熱段的中間，距加熱段底端115.6 cm的測量截面處于加熱段的出口位置，而距加熱段底端179.1 cm的測量截面處于加熱段下游的非加熱區(qū)。WARD 61棒束實(shí)驗(yàn)的具體幾何參數(shù)列于表3。圖6示出了WARD傾斜功率分布實(shí)驗(yàn)中加熱棒的功率在徑向上的功率分布。

圖5 WARD實(shí)驗(yàn)段組件子通道編號方案Fig.5 Sub-channel numbering scheme for WARD 61-rod bundle

表3 WARD 61棒束實(shí)驗(yàn)組件幾何參數(shù)Table 3 Geometric parameter of WARD 61-rod bundle experimental assembly

圖6 WARD實(shí)驗(yàn)組件加熱棒功率傾斜徑向分布 Fig.6 Inclined radial power distribution of heating rod of WARD experimental assembly

圖7示出了在徑向均勻功率分布下第243、218、227組實(shí)驗(yàn)不同軸向高度測量得到的無量綱溫度分布的計(jì)算結(jié)果與實(shí)驗(yàn)值。對于這3組實(shí)驗(yàn)，其功率與冷卻劑的流量大致呈正比，出口與入口的平均溫差約維持在100 ℃。表4示出了這3組實(shí)驗(yàn)的具體工況參數(shù)。在使用子通道分析方法來計(jì)算這3個工況時，湍流交混關(guān)系式選擇Zhukov模型。由于這3組實(shí)驗(yàn)的雷諾數(shù)范圍位于層流區(qū)與層流湍流過渡區(qū)，因此在阻力模型上選擇了Ki(Kirikkov)模型[26]與CT(Cheng-Todreas)[16]模型，這兩種模型均考慮了層流區(qū)的摩擦阻力系數(shù)以及層流湍流過渡區(qū)的插值。從圖7可看出，本文開發(fā)的子通道分析程序的計(jì)算結(jié)果與這幾組實(shí)驗(yàn)得出的結(jié)果符合得較好，并且選擇不同阻力模型帶來的不確定性較小。當(dāng)冷卻劑質(zhì)量流量減小時，不同的通道在橫向上得到的溫度會變得更加平滑。對于高流量時得到的溫度分布會有一個比較明顯的“拱形”形狀。這是由于邊通道的質(zhì)量流量高，因此邊通道的溫度相對較低。而通道間通過繞絲的交混以及分子-湍流擴(kuò)散傳遞能量，最終得到拱形的溫度分布。對于軸向高度為179.1 cm的測量截面，由于其位于非加熱區(qū)中，因此在非加熱區(qū)中通道間的能量傳遞過程使得通道間的溫差逐漸降低，其溫度分布相對于加熱段出口115.6 cm處的溫度分布更加平坦。對于低流量工況，通道間的能量交換過程主要依靠通道間的導(dǎo)熱過程。從實(shí)驗(yàn)227組的數(shù)據(jù)中可看出，在低流量工況下，179.1 cm處測量截面上的溫度分布幾乎是均勻的，在非加熱區(qū)中通道間的分子-湍流熱傳導(dǎo)效應(yīng)很快將溫度展平。

圖7 WARD均勻功率分布實(shí)驗(yàn)上、中、下測點(diǎn)的無量綱溫度分布計(jì)算值與實(shí)驗(yàn)值對比Fig.7 Comparison between calculated and experimental values of dimensionless temperature distribution of upper,middle and lower measurement points in WARD uniform power distribution experiment

表4 WARD 61棒束實(shí)驗(yàn)工況參數(shù)Table 4 WARD 61-rod bundle experiment operating condition parameter

圖8 WARD傾斜功率分布實(shí)驗(yàn)上、中、下測點(diǎn)的無量綱溫度分布計(jì)算值與實(shí)驗(yàn)值對比Fig.8 Comparison between calculated and experimental values of dimensionless temperature distribution of upper,middle and lower measurement points in WARD gradient power distribution experiment

圖8中同時示出了使用Cheng-Todreas繞絲棒束交混模型與使用Zhukov交混模型計(jì)算得到的各組實(shí)驗(yàn)的對比。從圖中可看出，相比于Zhukov交混模型，Cheng-Todreas交混模型在層流湍流過渡區(qū)223、221組實(shí)驗(yàn)得到的結(jié)果與實(shí)驗(yàn)值符合得更好。并且，在高功率分布一側(cè)的通道溫度相比Zhukov交混模型計(jì)算得到的值更低。由于冷卻劑在外圍間隙中的單向流動，較熱一側(cè)的邊通道的冷卻劑流向功率較低的一側(cè)，而較冷一側(cè)的邊通道的冷卻劑流向功率較高的一側(cè)。因此模型給出的外圍間隙單向流動強(qiáng)度越強(qiáng)，高功率分布一側(cè)的邊通道的溫度越低，而低功率分布的一側(cè)的邊通道溫度越高。由于Cheng-Todreas模型給出的繞絲帶來的外圍間隙的單向?qū)α鹘换鞆?qiáng)度要高于Zhukov模型，因此Cheng-Todreas模型計(jì)算得到的高功率側(cè)的邊通道99溫度更低。在層流湍流過渡區(qū)，Cheng-Todreas模型給出的通道間的交混似乎相比Zhukov能與實(shí)驗(yàn)符合得更好。但對于高流量工況313組，Cheng-Todreas模型預(yù)測得到的外圍間隙的單向橫流強(qiáng)度相比于真實(shí)值稍高，以至于高功率側(cè)的邊通道99計(jì)算得到溫度與實(shí)驗(yàn)值相比更低且低功率側(cè)的邊通道115得到的溫度相比實(shí)驗(yàn)值更高。對于加熱段中間的測量截面(圖中的藍(lán)線)，這兩種模型給出的溫度差別并不明顯。這是由于在距入口較近的軸向高度較低的位置，高低功率側(cè)的邊通道溫度相差并不太明顯，因此外圍間隙的強(qiáng)迫橫流強(qiáng)度對溫度分布影響不大。對于低流量實(shí)驗(yàn)229組，由于其溫度分布主要受通道簡單的熱傳導(dǎo)系數(shù)影響，因此兩種模型計(jì)算得到的溫度分布相差并不大。在總體上來說，Cheng-Todreas給出的溫度分布相比Zhukov模型與實(shí)驗(yàn)值符合得更好。

4 結(jié)論

本文利用雙區(qū)域模型開發(fā)了一款用于鈉冷快堆組件熱工水力分析的子通道分析程序SPLICA，SPLICA程序中通過雙區(qū)域模型本構(gòu)關(guān)系式詳細(xì)地考慮了繞絲帶來的通道間的交混效應(yīng)。通過與FFM-2A 19棒束實(shí)驗(yàn)以及WARD 61棒束實(shí)驗(yàn)溫度分布數(shù)據(jù)進(jìn)行對比，驗(yàn)證了本文開發(fā)的子通道程序在層流、湍流以及層流湍流過渡區(qū)工況下對鈉冷快堆組件熱工水力分析具有良好的適用性，并且在組件徑向功率傾斜情況下計(jì)算結(jié)果也具有較高的準(zhǔn)確度。對于FFM 2A實(shí)驗(yàn)，與經(jīng)過二次開發(fā)的COBRA-Ⅳ程序相比，本文開發(fā)的子通道分析程序與實(shí)驗(yàn)值符合得更好。對于WARD 61棒束傾斜功率分布實(shí)驗(yàn)，Cheng-Todreas交混模型在層流湍流過渡區(qū)計(jì)算得到的溫度分布較Zhukov交混模型與實(shí)驗(yàn)符合得更好。本程序能為池式鈉冷快堆組件的熱工水力研究提供有效的設(shè)計(jì)和分析工具。