師 泰，張東輝，劉一哲

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

示范快堆是我國自主設(shè)計建造的一座池式鈉冷快中子反應(yīng)堆，反應(yīng)堆采用六角形組件作為燃料，燃料在反應(yīng)堆運行過程中浸泡在導(dǎo)熱性能較好的液態(tài)金屬鈉中，在換料過程中會通過換熱性能較差的氬氣空間。因此，保證燃料組件在氬氣環(huán)境中安全性是反應(yīng)堆安全設(shè)計的重要組成部分[1]?？於讶剂辖M件初始釋熱較大，一般采用乏燃料組件裝入堆內(nèi)儲存阱中放置1～2個換料周期的方式降低其衰變熱功率。換料開始后，組件由換料系統(tǒng)提升到反應(yīng)堆外部的轉(zhuǎn)運室中，轉(zhuǎn)運室為氬氣環(huán)境，長時間懸停在轉(zhuǎn)運室可能導(dǎo)致組件超溫，甚至破損，組件破損熔化會導(dǎo)致放射性物質(zhì)進入環(huán)境[2-3]。因此，研究燃料組件在氬氣環(huán)境下的換熱特性，保證反應(yīng)堆換料工程中組件的安全極其重要。

目前國內(nèi)外針對六角形組件在氣體環(huán)境下的換熱特性研究甚少，美國桑迪亞國家實驗室[4]開展了模擬水平運輸過程快堆217棒乏燃料組件的換熱特性試驗，進行了低加熱功率范圍在不同氣體介質(zhì)中的傳熱試驗，得到了組件內(nèi)部溫度分布并擬合了熱導(dǎo)率經(jīng)驗關(guān)系式[5]。Alyokhina等[4-6]建立含多個快堆乏燃料組件的儲存桶的三維模型，研究環(huán)境溫度變化對貯存桶內(nèi)溫度分布的影響，結(jié)果表明外部環(huán)境溫度變化對貯存桶內(nèi)溫度分布無明顯影響，同時解決了反向共軛傳熱問題。

本文采用三維CFD軟件研究37棒乏燃料組件在充滿氬氣環(huán)境的轉(zhuǎn)運室懸停時，燃料組件依靠自然循環(huán)方式冷卻的穩(wěn)態(tài)溫度場，同時采用37棒模擬組件實驗進行驗證。

1 實驗簡介

該實驗?zāi)M37棒組件在氬氣環(huán)境下的換熱，試驗裝置如圖1所示。試驗段懸掛在壓力容器中，通過對試驗段組件進行電加熱模擬組件的發(fā)熱，組件內(nèi)布置有熱電偶，通過熱電偶測量組件內(nèi)的溫度分布，從而得到組件在不同功率下的局部溫度分布。壓力容器采用三段式設(shè)計，氬氣容器采用真空泵抽真空后注入氬氣，氬氣壓力穩(wěn)定在(0.10±0.01) MPa，為滿足容器壁面穩(wěn)定的溫度邊界條件，在壓力容器外壁焊接有控制邊界溫度冷卻管路。壓力容器內(nèi)徑為500 mm，分3段設(shè)計，下段高度為400 mm，中段高度為1 200 mm，上段高度為800 mm。

圖1 試驗裝置示意圖Fig.1 Schematic of experimental apparatus

圖2 試驗段Fig.2 Test section

37棒模擬組件試驗段豎直安裝于壓力容器中段中心位置，試驗段由棒束和組件盒組成，如圖2所示。棒束呈六角形排列，由30根加熱棒和7根測溫棒組成。加熱棒的內(nèi)部結(jié)構(gòu)如圖3所示，棒內(nèi)使用導(dǎo)熱性能較好的MgO粉末填充壓實。根據(jù)導(dǎo)電率隨溫度的變化關(guān)系，將處于對稱位置的加熱棒組合為5組，即H1～H5，如圖2b所示。每組加熱棒并聯(lián)后與程控電源連接，獨立調(diào)節(jié)電源電壓，保證加熱功率均勻。

組件盒為正六邊形結(jié)構(gòu)，內(nèi)對邊距為44 mm，厚度為3 mm。元件棒直徑為6 mm，棒中心距為6.95 mm，總長度為1.2 m。

圖3 加熱棒的結(jié)構(gòu)Fig.3 Structure of heated rod

中心測溫棒(M1)內(nèi)安裝4只熱電偶，外圈測溫棒(M2～M7)分別安裝2只熱電偶，測溫棒不具備加熱功能，測溫棒的測點位置如圖4所示。

圖4 M1～M7內(nèi)熱電偶位置Fig.4 Thermocouple locations in M1-M7

2 計算程序及模型

2.1 計算程序

本研究采用 CD-asapco公司開發(fā)的STAR-CCM+軟件進行計算，該軟件采用最先進的連續(xù)介質(zhì)力學(xué)數(shù)值技術(shù)開發(fā)的新一代CFD求解器。它搭載了CD-asapco獨創(chuàng)的最新網(wǎng)格生成技術(shù)，可完成復(fù)雜形狀數(shù)據(jù)的輸入。它的主要特點在于條件設(shè)置與后處理方面，如同時打開兩個算例，某些條件設(shè)置僅需復(fù)制粘貼即可互相轉(zhuǎn)換，省去了再輸入1遍的時間，尤其是在氣體吸收系數(shù)各波段的波長范圍錄入以及吸收系數(shù)錄入等數(shù)據(jù)量大的設(shè)置中，為用戶節(jié)省了大量時間。另一方面，STAR-CCM+在多面體網(wǎng)格生成方面有較大的優(yōu)勢，尤其是應(yīng)對復(fù)雜幾何結(jié)構(gòu)的計算，如燃料棒束、堆芯等復(fù)雜幾何體[7-9]。

2.2 關(guān)鍵模型分析

1) 自然循環(huán)能力估算

熱能傳遞的方式有3種：熱傳導(dǎo)、熱對流和熱輻射。乏燃料在氬氣環(huán)境下沒有強制對流，組件內(nèi)部靠自然循環(huán)建立流動特性，由于組件內(nèi)部間隙較小，自然循環(huán)能力低，可忽略不計。自然循環(huán)能力估算如下。

假設(shè)組件與氬氣充分換熱，氬氣出口溫度與組件溫度一致，元件棒表面溫度為700 ℃，氬氣在組件內(nèi)的自然循環(huán)流量計算公式[10-11]為：

(1)

其中：qm為流體質(zhì)量流量，kg/s；ρ為密度，kg/m3；g為重力加速度，m/s2；A為流體的橫截面積，m2；αv為體膨脹系數(shù)，K-1；tw為組件表面溫度，℃；ta為氬氣環(huán)境溫度(與環(huán)境溫度一致，50 ℃)；δ為元件棒之間的間隙；ν為氣體的運動黏度，m2/s。

氬氣自然循環(huán)所帶走的熱量為：

Q=Δt·qm·cp

(2)

其中：Q為熱量，W；Δt為氬氣溫度變化量，℃；cp為比定壓熱容，J/(kg·℃)。計算得到Q為17.187 W，相對于400～1 000 W組件功率可忽略不計。

2) 繞絲的模型簡化

由于組件內(nèi)部流動特性較小，繞絲對組件換熱特性的影響主要由于繞絲較好的導(dǎo)熱性能，而繞絲對組件內(nèi)部流動的影響可忽略。因此，CFD建模中采用無繞絲的模型，繞絲對組件內(nèi)部導(dǎo)熱的影響采用等效導(dǎo)熱的方法分析，無繞絲的等效導(dǎo)熱方法相對于真實組件的計算結(jié)果保守。采用等效熱阻方法得到等效導(dǎo)熱系數(shù)[12-14]如下：

(3)

其中：λt為等效導(dǎo)熱系數(shù)；λa和λs分別為氬氣和鋼的導(dǎo)熱系數(shù),W/(m·K)；Aa和As分別為氬氣和鋼的等效面積,m2。氬氣的導(dǎo)熱系數(shù)為0.036 W/(m·K)，計算得到等效導(dǎo)熱系數(shù)為0.051 W/(m·K)。

3) 輻射模型

S2S(surface-to-surface)輻射模型適用于計算封閉系統(tǒng)中各反射面之間的輻射傳熱，并且模型主要參數(shù)運行在計算迭代步數(shù)開始前預(yù)處理完成，因此計算具有較高的效率，故本文選擇S2S模型[15]。

對空間兩個微元面而言，其輻射換熱如圖5所示，由dS1面發(fā)射被dS2面吸收的輻射總功率為：

P1-2=i′1dS1cosβ1(dS2cosβ2/L2)

(4)

其中：P為功率；i′為輻射強度，W/(m2·sr)；S為面積；L為特征長度；β為空間角度。

圖5 微元面輻射換熱Fig.5 Micro-element radiative heat transfer

兩微元面之間的輻射角系數(shù)可由下式計算獲得：

(5)

其中，F(xiàn)為角系數(shù)。當兩個黑體表面在相同溫度場中達到熱平衡時，其相應(yīng)的輻射角系數(shù)應(yīng)滿足如下關(guān)系式：

dFi-jdSi=dFj-idSj

(6)

據(jù)此拓展到空間有限大兩平面之間時，有：

(7)

同理可知：

(8)

2.3 邊界條件設(shè)置及計算工況

乏燃料組件在轉(zhuǎn)運室的換熱是一個自然對流冷卻的過程，由于自然循環(huán)能力較弱，因此忽略對流的影響，組件主要通過導(dǎo)熱和輻射散熱，計算為穩(wěn)態(tài)計算。邊界條件主要包括熱源邊界、固體邊界和外壁面邊界。

1) 熱源邊界，乏燃料組件在出堆后組件功率軸向分布可近似為恒定值，因此元件棒功率采用體積功率的輸入?yún)?shù)，對30根元件棒軸向0.2～1 m段賦予相同的體積功率。實驗中加熱棒為30根元件棒，因此400、800和1 000 W功率下體功率分別為5.89×105、1.18×106、1.47×106W/m3。

2) 固體邊界，乏燃料組件內(nèi)部元件棒與氬氣和氬氣與組件盒之間為接觸界面邊界條件。固體導(dǎo)熱系數(shù)依據(jù)指定的導(dǎo)熱系數(shù)進行計算，固體之間無接觸熱阻。輻射傳輸采用氬氣內(nèi)部輻射傳輸。

3) 外壁面邊界，外壁面邊界條件主要包括組件盒外壁和上下端面。該邊界條件采用壁面邊界條件，由于實驗過程中組件出現(xiàn)部分氧化現(xiàn)象，組件盒外壁的發(fā)射率如圖6所示，表面對流換熱系數(shù)均采用5 W/(m2·K)，環(huán)境溫度為常數(shù)50 ℃。上下端面對組件換熱影響較小，故采用絕熱壁面邊界條件。

圖6 歸一化發(fā)射率分布Fig.6 Normalized emissivity distribution

CFD計算與實驗的工況相同，計算輸入?yún)?shù)為：氬氣導(dǎo)熱系數(shù)，0.051 W/(m·K)；組件盒和元件棒導(dǎo)熱系數(shù)，16 W/(m·K)；表面換熱系數(shù)，5 W/(m2·K)；環(huán)境溫度，50 ℃；功率，400、800、1 000 W。氣體的等效導(dǎo)熱率為0.051 W/(m·K)。

2.4 網(wǎng)格劃分及收斂性

多面體網(wǎng)格相對于四面體網(wǎng)格能在復(fù)雜幾何體及邊界層、長通道或小間隙的區(qū)域可在較少的控制體的基礎(chǔ)上達到較高的精度。同時，由于多面體網(wǎng)格有很多鄰居單元，所以能精確計算控制體的梯度。但多面體網(wǎng)格可能造成較大的計算量和內(nèi)存需求。由于組件尺寸跨度大、結(jié)構(gòu)復(fù)雜，采用多面體網(wǎng)格可提高計算精度。整體網(wǎng)格如圖7所示，網(wǎng)格質(zhì)量列于表1。表1中，網(wǎng)格在拓撲上有效，并且沒有負體積，總網(wǎng)格數(shù)為5 035 336。

圖7 37棒組件網(wǎng)格圖Fig.7 37 rods assembly mesh

表1 網(wǎng)格質(zhì)量Table 1 Mesh quality

本文分別采用200萬、400萬、500萬、900萬和2 200萬網(wǎng)格進行了測試計算，中心棒600 mm位置溫度如圖8所示。結(jié)果表明，500萬網(wǎng)格與2 200萬網(wǎng)格對應(yīng)的組件中心棒600 mm位置溫度的相對誤差為0.29%，能滿足要求，故采用500萬網(wǎng)格進行分析。

由于乏燃料組件僅考慮導(dǎo)熱和輻射，組件盒內(nèi)氬氣采用固體模型，計算收斂性主要參考能量殘差和局部溫度恒定。圖9示出計算迭代過程中的殘差和中心棒800 mm處的溫度隨迭代步數(shù)的變化，圖中可看到迭代步數(shù)4 500步之后組件溫度基本不變，能量殘差小于10-5，故判斷計算收斂。

圖8 網(wǎng)格數(shù)與測量值的關(guān)系Fig.8 Relationship between number of grid and measured value

圖9 能量殘差和溫度隨迭代步數(shù)的變化Fig.9 Variation of energy residual and temperature with number of iteration step

3 實驗和計算結(jié)果分析

乏燃料組件溫度直接關(guān)系到乏燃料組件操作的設(shè)計及安全性，尤其是燃料元件包殼的最高溫度，是關(guān)系到組件能否包容放射性的重要因素。數(shù)值模擬結(jié)果與實驗測量值相對誤差較小，能滿足CFD模擬實驗的精度要求。由于采用的填充材料為氧化鎂，該材料導(dǎo)熱性能較好，實驗中，中心棒600～800 mm的溫差較小(<6 ℃)，可近似認為實驗和計算的中心棒600、650、800 mm 3個位置的最高溫度為組件最高溫度。

實驗和數(shù)值模擬結(jié)果表明中心測溫棒的溫度最高，且最高溫度出現(xiàn)在600～800 mm之間，400 W功率下最高溫度為440 ℃，800 W功率下最高溫度為581 ℃，1 000 W功率下最高溫度為631 ℃。實驗和數(shù)值模擬不同功率下中心棒M1最高溫度對比如圖10a所示，實驗和數(shù)值模擬不同功率下盒壁最高溫度對比如圖10b所示，實驗和數(shù)值模擬中心棒M1的溫度分布對比如圖11所示。

圖10 不同功率下中心棒M1和盒壁最高溫度對比Fig.10 Comparison of the highest temperatures for central rod M1 and component box at different powers

圖11 不同功率下中心棒M1溫度對比Fig.11 Comparison of temperature for central rod M1 at different powers

不同功率下中心棒M1溫度CFD計算值和實驗值相對誤差對比如圖12所示。由圖12可知，實驗數(shù)據(jù)和CFD模擬的相對誤差較小，最大的相對誤差為6%，出現(xiàn)在距離上下端較近位置，這是由于在實驗中采用鋼條對實驗裝置進行了固定，鋼條較好的導(dǎo)熱性能對實驗產(chǎn)生較大的影響。CFD計算中最高溫度能較為精確地反映實驗中最高溫度，中心棒600 mm和650 mm位置的溫度相對誤差不超過1%。

圖12 不同功率下中心棒M1溫度的相對誤差Fig.12 Relative error of temperature for central rod M1 at different powers

實驗和數(shù)值模擬結(jié)果表明外圈組件溫度較中心棒組件溫度低，由于實驗和CFD計算中外圈測溫棒不同位置的測溫點不在同一元件棒上，由于不同元件棒相對于中心棒的角系數(shù)有一定的區(qū)別，輻射換熱傳遞的熱量不同導(dǎo)致溫度出現(xiàn)較小的波動。外圈組件的最高溫度出現(xiàn)在600～800 mm之間，CFD模擬中400 W功率下最高溫度為407 ℃，800 W功率下最高溫度為537 ℃，1 000 W功率下最高溫度為587 ℃。實驗和數(shù)值模擬外圈測溫棒M2的溫度分布對比如圖13所示。不同功率下外圍組件溫度CFD計算值和實驗值相對誤差對比如圖14所示。由圖14可知，實驗數(shù)據(jù)和CFD模擬的相對誤差較小，最大的相對誤差為3%。

圖13 不同功率下外圈測溫棒M2溫度對比Fig.13 Comparison of temperature for outer test rod M2 at different powers

圖14 不同功率下外圍測溫棒M2溫度的相對誤差Fig.14 Relative error of temperature for outer test rod M2 at different powers

實驗和數(shù)值模擬結(jié)果表明盒壁溫度較燃料元件溫度低，由于實驗和CFD中盒壁測溫點處于六角形各邊的中心，測量溫度為盒壁的最高溫度，盒壁溫度在6個角位置較低。盒壁最高溫度出現(xiàn)在600～800 mm之間，CFD模擬中400 W功率下最高溫度為285 ℃，800 W功率下最高溫度為373 ℃，1 000 W功率下最高溫度為405 ℃。實驗和CFD模擬盒壁溫度分布對比如圖15所示。不同功率下外圍組件溫度CFD計算值和實驗值相對誤差對比如圖16所示。由圖16可知，實驗數(shù)據(jù)和CFD模擬的相對誤差較小，最大的相對誤差為9%。

圖15 不同功率下盒壁溫度對比Fig.15 Comparison of temperature for component box at different powers

圖16 不同功率下盒壁溫度的相對誤差Fig.16 Relative error of temperature for component box at different powers

4 結(jié)論

本文采用實驗和CFD模擬兩種方法對37棒乏燃料組件在氬氣環(huán)境下的換熱特性進行了研究，研究結(jié)果表明：

1) 數(shù)值模擬采用的方法能較好模擬實驗工況，中心棒溫度最大相對誤差為6%，最高溫度的相對誤差小于1%，實驗和數(shù)值模擬結(jié)果表明組件峰值溫度出現(xiàn)在中心棒600～800 mm之間。

2) 乏燃料在氬氣環(huán)境下的自然循環(huán)能力較小，乏燃料主要通過氬氣的熱傳導(dǎo)和輻射向外界傳遞熱量，其中組件在400、800和1 000 W功率下輻射換熱占總換熱量的比例分別為37%、52%和56%。

3) 37棒六角形的乏燃料組件在0.5、1、1.25 kW/m的線功率密度下，穩(wěn)態(tài)計算最高溫度分別為440、581、631 ℃，能保證組件的安全性。