明 春，韓智杰，何曉軍，初泉麗

(1.中國(guó)原子能科學(xué)研究院 反應(yīng)堆工程技術(shù)研究部，北京 102413；2.國(guó)家核安保技術(shù)中心，北京 102401)

利用燃料組件導(dǎo)向管的內(nèi)部空間開(kāi)展材料或釋熱元件的輻照考驗(yàn)，是研究材料輻照性能及釋熱元件安全特性較常用的一種輻照試驗(yàn)方法。以典型壓水堆燃料組件為例，燃料棒按照正方形柵格的形式排列，外側(cè)為8根普通燃料棒，中心為導(dǎo)向管，內(nèi)部插入釋熱元件棒。該結(jié)構(gòu)中存在棒束通道外流流道-導(dǎo)向管-窄環(huán)縫內(nèi)流流道-釋熱元件的復(fù)雜傳熱通道，其中，導(dǎo)向管內(nèi)熱流密度Q0未知，且內(nèi)外流流道差異較大，實(shí)際計(jì)算時(shí)涉及到多個(gè)熱邊界的溫度耦合問(wèn)題，需聯(lián)立求解多組傳熱微分方程，現(xiàn)有程序通常不會(huì)單獨(dú)考慮這一特殊結(jié)構(gòu)，難以計(jì)算出精確結(jié)果。工程上一般采用CFD軟件進(jìn)行仿真模擬，需進(jìn)行三維建模，耗時(shí)較長(zhǎng)、流程繁瑣。因此，本文擬設(shè)計(jì)專(zhuān)用的耦合算法用于解決套管結(jié)構(gòu)的傳熱問(wèn)題，以便直接通過(guò)程序計(jì)算出符合燃料性能分析程序需要的套管結(jié)構(gòu)溫度數(shù)據(jù)。

實(shí)際計(jì)算時(shí)需要相應(yīng)的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行傳熱關(guān)系式擬合。Weisman[1]和Markoczy[2]分別針對(duì)全堆芯和堆內(nèi)有限區(qū)域給出了相應(yīng)的傳熱關(guān)系式。對(duì)于環(huán)管流道，孫立成等[3]、曾和義等[4]和白博峰等[5]均開(kāi)展過(guò)環(huán)管加熱實(shí)驗(yàn)，得到了不同加熱條件下環(huán)管傳熱的擬合關(guān)系式，Liu等[6]對(duì)環(huán)管傳熱的離心效應(yīng)進(jìn)行了深入研究。

本文擬基于以上文獻(xiàn)，依據(jù)流場(chǎng)基本情況，通過(guò)迭代算法和熱工水力基本理論，開(kāi)發(fā)具有外套管的釋熱元件表面溫度分布計(jì)算程序，并與燃料性能分析程序進(jìn)行耦合，為具有外套管的釋熱元件的設(shè)計(jì)及安全評(píng)價(jià)提供技術(shù)支持。

1 幾何模型

根據(jù)對(duì)稱(chēng)性，選擇典型求解單元(圖1虛線內(nèi)范圍)作為研究對(duì)象，不考慮組件所處具體堆芯位置對(duì)溫度計(jì)算的影響。

導(dǎo)向管外側(cè)流體所處流道為典型的棒束流道，本文將外流流道等效處理為圓形流道，則整個(gè)求解單元可視為1個(gè)多層套管結(jié)構(gòu)，如圖2所示，最外側(cè)為燃料棒外壁，向內(nèi)分別為外流流道、導(dǎo)向管、內(nèi)流流道和釋熱元件。

為與性能分析程序計(jì)算節(jié)點(diǎn)相匹配，網(wǎng)格劃分遵循1.5維的處理方式，首先將多層套管求解單元在軸向劃分為若干個(gè)短圓柱段，對(duì)于每個(gè)軸向圓柱段，徑向存在不同的傳熱部件。具體節(jié)點(diǎn)劃分示于圖3。

圖1 求解單元結(jié)構(gòu)示意圖Fig.1 Structure of solving unit

圖2 多層套管結(jié)構(gòu)示意圖Fig.2 Structure of multi-casing

圖3 求解域節(jié)點(diǎn)劃分Fig.3 Meshing of solution region

2 理論模型

2.1 冷卻劑軸向傳熱

對(duì)于流體區(qū)域，已知兩側(cè)邊界的熱流密度，根據(jù)能量守恒關(guān)系有：

qΔl=cpqm(To-Ti)

(1)

其中：q為流體邊界上線熱流密度的代數(shù)和；Δl為換熱邊界長(zhǎng)度，即軸向段長(zhǎng)度；cp為比定壓熱容；qm為質(zhì)量流量；To和Ti分別為軸向段出、入口溫度，定性溫度取軸向段流體平均溫度。

2.2 外流流道對(duì)流換熱

導(dǎo)向管外側(cè)流體所處流道為典型的棒束流道[6]，在進(jìn)行單向流分析時(shí)，需單獨(dú)考慮通道形狀對(duì)傳熱系數(shù)的影響。實(shí)際工程計(jì)算中，一般先采用Dittus-Boelter關(guān)系式[7]計(jì)算得到等效圓管的努塞爾數(shù)Nu∞,cir，然后進(jìn)行系數(shù)修正得到棒束通道的Nu，即：

Nu=ψNu∞,cir

(2)

其中，ψ為修正系數(shù)。

對(duì)于正方形排列的棒束，考慮采用文獻(xiàn)[8]給出的通用關(guān)系式計(jì)算ψ：

ψ=1+0.912Re-0.1Pr0.4(1-2.004 3e-B)

(3)

(4)

2.3 內(nèi)流對(duì)流換熱

估算可知，內(nèi)流通道內(nèi)Re<10 000，而Dittus-Boelter關(guān)系式僅適用于旺盛湍流換熱，因此Nu需采用適用于過(guò)渡區(qū)的Gnielinski關(guān)系式[9]計(jì)算：

(5)

其中，下標(biāo)f和w分別表示以流體溫度和壁面溫度作為定性溫度。

根據(jù)白博峰等[5]的實(shí)驗(yàn)，環(huán)形管內(nèi)的湍流換熱相對(duì)于普通圓管內(nèi)湍流換熱有所加強(qiáng)。本文采用白博峰等進(jìn)行的單側(cè)加熱環(huán)管實(shí)驗(yàn)得到的擬合關(guān)聯(lián)式作為內(nèi)流Nu的計(jì)算式：

Nu=0.023Re0.91Pr0.4

(6)

式(6)的適用范圍為2 300

2.4 導(dǎo)向管導(dǎo)熱

對(duì)于導(dǎo)向管區(qū)域，可視為無(wú)內(nèi)熱源的一維環(huán)形構(gòu)件導(dǎo)熱問(wèn)題，導(dǎo)向管穩(wěn)態(tài)傳熱方程[10]為：

(7)

其中：r為節(jié)點(diǎn)半徑；T為節(jié)點(diǎn)溫度；kc為導(dǎo)向管熱導(dǎo)率。定性溫度取導(dǎo)向管平均溫度。

3 計(jì)算方法

穩(wěn)態(tài)情況下，傳熱區(qū)域的邊界上熱流密度恒定，即導(dǎo)向管兩側(cè)熱流密度Q0相同。通過(guò)求解Q0，可使內(nèi)外流流道傳熱邊界條件封閉，從而完成傳熱方程的求解。

本文所設(shè)計(jì)算法以軸向段作為基本求解單元，將Q0作為收斂標(biāo)準(zhǔn)，首先通過(guò)初始Q0以及傳熱相關(guān)方程逐步求解出外流流道、導(dǎo)向管和內(nèi)流流道相應(yīng)的節(jié)點(diǎn)溫度，然后通過(guò)能量守恒反推出新的Q0，與初始值進(jìn)行收斂判斷，不收斂則重復(fù)上述過(guò)程直至收斂。Q0收斂后利用對(duì)流換熱關(guān)系可求解出釋熱元件最外側(cè)節(jié)點(diǎn)溫度，完成當(dāng)前軸向段溫度計(jì)算，并將結(jié)果保存，作為下一軸向段計(jì)算的邊值條件，最終完成全長(zhǎng)溫度計(jì)算。具體計(jì)算流程示于圖4。

圖4 軸向段溫度求解流程Fig.4 Solution procedure of axial segment temperature

本文采用Fortran語(yǔ)言，依據(jù)此算法編譯完成套管結(jié)構(gòu)溫度計(jì)算程序，用于對(duì)具有外套管的釋熱元件進(jìn)行全長(zhǎng)溫度計(jì)算。

4 計(jì)算結(jié)果及對(duì)比

利用根據(jù)上述算法所建立的套管結(jié)構(gòu)溫度分布計(jì)算程序，采用典型壓水堆燃料輸入?yún)?shù)(表1)和ANSYS Fluent軟件對(duì)程序進(jìn)行對(duì)比分析。

表1 程序輸入?yún)?shù)Table 1 Input parameters of program

4.1 計(jì)算結(jié)果

采用本文所建程序計(jì)算的套管結(jié)構(gòu)內(nèi)各區(qū)域的節(jié)點(diǎn)溫度示于圖5。由圖5可見(jiàn)，外流流道與導(dǎo)向管外壁溫度隨高度的增加而上升，而導(dǎo)向管內(nèi)壁、內(nèi)流流道、釋熱元件外壁溫度隨高度的增加呈先上升后緩慢下降的趨勢(shì)，溫度峰值出現(xiàn)在出口附近。

圖5 溫度分布計(jì)算結(jié)果Fig.5 Calculation result of temperature distribution

4.2 結(jié)果對(duì)比

采用ANSYS Fluent軟件對(duì)圖1的3×3棒束結(jié)構(gòu)建立三維模型，采用表1的邊界參數(shù)進(jìn)行熱工水力仿真計(jì)算，并與本文程序計(jì)算結(jié)果進(jìn)行對(duì)比，重點(diǎn)關(guān)注內(nèi)、外流流道和釋熱元件包殼外側(cè)的溫度分布情況，以及導(dǎo)向管兩側(cè)的熱流密度，結(jié)果示于圖6、7。由圖6、7可見(jiàn)，二者計(jì)算結(jié)果偏差較小，F(xiàn)luent仿真模擬所得溫度略高于程序計(jì)算值，相對(duì)誤差在5%以內(nèi)，本程序可在工程設(shè)計(jì)中提供一定參考。

圖6 Q0本文程序計(jì)算值和Fluent模擬值的對(duì)比Fig.6 Comparison of Q0 calculated by code in this paper and simulated from Fluent

5 程序應(yīng)用

套管結(jié)構(gòu)溫度計(jì)算程序改寫(xiě)為單獨(dú)的計(jì)算模塊后，可與傳統(tǒng)性能分析程序相耦合，在性能分析計(jì)算開(kāi)始前，調(diào)用該計(jì)算模塊完成釋熱元件表面溫度計(jì)算，結(jié)果作為邊界條件在后續(xù)計(jì)算中使用。

圖7 溫度分布本文程序計(jì)算結(jié)果與Fluent模擬值的對(duì)比Fig.7 Comparison of temperature distribution calculated by code in this paper and simulated from Fluent

采用已耦合套管溫度計(jì)算模塊的燃料性能分析程序，進(jìn)行完整的全壽期燃料性能分析，實(shí)現(xiàn)對(duì)燃料溫度、變形、裂變氣體釋放等行為的研究[11-13]。通過(guò)假想算例，研究核電站全壽期內(nèi)燃料性能的變化情況，尺寸參數(shù)及燃料棒平均線功率與表1一致，但釋熱元件線功率會(huì)隨時(shí)間推移逐漸下降，因此其全壽期平均值設(shè)定為1 kW/m。

釋熱元件典型性能分析結(jié)果示于圖8～10。由圖8可見(jiàn)，壽期初芯塊溫度相對(duì)較高，溫度分布曲線較陡峭，壽期末溫度普遍降低，分布趨于平緩[14]。分析圖9可知，前期由于裂變氣體持續(xù)穩(wěn)定釋放，氣體內(nèi)壓增大，而后期元件功率下降，裂變氣體釋放受到影響[15]，內(nèi)壓增長(zhǎng)略有減緩。由圖10可知，包殼與芯塊的徑向膨脹基本穩(wěn)定，軸向伸長(zhǎng)趨勢(shì)則逐步減緩。

圖8 壽期初與壽期末芯塊軸向溫度分布Fig.8 Axial temperature distribution of pellet during BOL and EOL

圖9 壽期內(nèi)元件內(nèi)壓的變化Fig.9 Variation of element inner pressure during lifetime

圖10 壽期內(nèi)元件力學(xué)變形的變化Fig.10 Variation of element mechanical deformation during lifetime

綜上所述，基于本文算法修改的燃料性能分析程序在完整實(shí)現(xiàn)原有性能分析程序計(jì)算功能的前提下，可自行計(jì)算燃料棒外溫度分布情況，并將傳熱方程求解所需邊界條件傳遞給溫度模塊，保證后續(xù)運(yùn)算流程正常進(jìn)行，從而避免了對(duì)熱工水力程序的依賴(lài)，提高了分析流程的獨(dú)立性。

6 結(jié)論

本文針對(duì)套管傳熱結(jié)構(gòu)設(shè)計(jì)了專(zhuān)門(mén)的迭代算法用于域內(nèi)溫度分布求解，并基于該算法開(kāi)發(fā)了相應(yīng)的Fortran程序模塊，實(shí)現(xiàn)了具有外套管的釋熱元件包殼溫度的求解。與大型CFD商用軟件Fluent的仿真模擬結(jié)果比較，兩者差異在5%以內(nèi)，取得了較好的一致性。

耦合套管結(jié)構(gòu)傳熱計(jì)算模塊的燃料性能分析程序可自行完成溫度邊界條件的計(jì)算，獨(dú)立完整地進(jìn)行性能分析流程，減少了計(jì)算流程的復(fù)雜度。目前，該程序已實(shí)現(xiàn)工程應(yīng)用，相關(guān)釋熱元件已入堆開(kāi)展輻照試驗(yàn)。