基于SANA實(shí)驗(yàn)的三維熱工分析程序計(jì)算驗(yàn)證

張向丞, 佘頂, 陳福冰, 石磊

(清華大學(xué) 核能與新能源技術(shù)研究院, 北京 100084)

球床式高溫氣冷堆是一種具有第4代核能系統(tǒng)特征的先進(jìn)反應(yīng)堆,其固有安全性能夠確保在各種事故工況下,堆芯的燃料最高溫度不會超過設(shè)計(jì)限值,從而避免大量放射性泄漏[1]。二維系統(tǒng)分析程序THERMIX和TINTE等被廣泛應(yīng)用于高溫氣冷堆內(nèi)的熱工分析,但這些程序的建模和計(jì)算都是基于二維的圓柱對稱假設(shè),難以處理高溫氣冷堆內(nèi)存在的由于功率分布不均勻或幾何結(jié)構(gòu)不對稱引起的三維熱工問題。因此,開展針對球床式高溫氣冷堆的三維熱工水力分析十分必要[2]。商用流體動力學(xué)計(jì)算軟件如FLUENT和CFX等可以應(yīng)用于球床高溫氣冷堆的局部三維熱工水力分析,但全堆建模的時效性、復(fù)雜性以及難以與中子動力學(xué)計(jì)算程序進(jìn)行耦合等問題仍然存在。因此,開發(fā)能夠?qū)崿F(xiàn)全堆快速三維熱工計(jì)算的系統(tǒng)分析程序?qū)τ谙乱淮邷貧饫涠训脑O(shè)計(jì)工作具有重要價(jià)值[3]。

二維熱工分析程序DAYU2D在吸收THERMIX軟件包的基礎(chǔ)上開發(fā)。基于該程序框架,本文根據(jù)堆芯內(nèi)部的守恒方程和適用于球床式高溫氣冷堆的相關(guān)經(jīng)驗(yàn)關(guān)系式,建立了三維熱工水力計(jì)算模型,開發(fā)了三維熱工分析程序DAYU3D,并針對德國SANA基準(zhǔn)實(shí)驗(yàn)的部分工況進(jìn)行了計(jì)算分析,驗(yàn)證了程序的計(jì)算功能。

1 DAYU3D程序三維熱工分析模型

DAYU2D程序是由清華大學(xué)核能與新能源技術(shù)研究院開發(fā)的二維熱工分析程序,用于分析球床式高溫氣冷堆熱工水力特性。該程序采用了THERMIX程序中驗(yàn)證過的模型,包括固體材料和冷卻劑的守恒方程以及材料的物性等。目前,DAYU2D程序已經(jīng)完成了HTR-10和HTR-PM的穩(wěn)態(tài)和瞬態(tài)算例驗(yàn)證,取得了良好的效果。在該程序框架下,本文開發(fā)的DAYU3D從基本守恒方程出發(fā),在三維柱坐標(biāo)系下對守恒方程進(jìn)行簡化、離散,結(jié)合球床傳熱及流動的經(jīng)驗(yàn)關(guān)系式,建立了三維熱工計(jì)算模型[4]。

1.1 三維熱傳導(dǎo)計(jì)算模型

球床中的固體區(qū)域能量守恒方程,考慮體積為dV的控制體,可以表示為:

(1)

式中:下標(biāo)s為固體;cv,s、Ts、qn、λeff,s、qk分別為固體的定容比熱、溫度、核熱源密度、球床等效導(dǎo)熱系數(shù)以及通過對流換熱進(jìn)入到控制體的熱量率,其中λeff,s由Zehner-Schlünder關(guān)系確定[5]。

圖1為固體節(jié)點(diǎn)(I,J,K)處的計(jì)算網(wǎng)格劃分,在對上述方程進(jìn)行離散時,空間上采用了中心差分格式,時間離散選用后向差分,差分網(wǎng)格的邊界則由周圍的材料網(wǎng)格的體積中心連接構(gòu)成,待求解的固體溫度位于材料網(wǎng)格的端點(diǎn)。

圖1 固體節(jié)點(diǎn)(I,J,K)處的差分網(wǎng)格Fig.1 Differential grid at solid node (I, J, K)

由于球床在進(jìn)行熱工分析時是被視作各向同性的多孔介質(zhì)模型,因此相較于二維模型,三維熱傳導(dǎo)計(jì)算模型只需多考慮周向上的熱量傳遞[6]。

1.2 三維對流計(jì)算模型

氣體相關(guān)參數(shù)的求解采用多孔介質(zhì)模型,并且假設(shè)氣體為理想流體,吸收的熱量均轉(zhuǎn)化為焓,滿足h=cp,gTg。可得關(guān)于流體的能量守恒方程:

(λeff,gΔTg)+αF(Ts-Tg)

(2)

式中:下標(biāo)g為氣體;ρg為氣體密度;cp,g為定壓比熱;Tg為氣體溫度;v為氣體的速度矢量;λeff,g為氣體有效導(dǎo)熱系數(shù);α、F分別為氣體和固體間的對流換熱系數(shù)以及對流換熱面積。

氣體的質(zhì)量守恒方程可以表示為:

(3)

式中G為氣體的質(zhì)量流密度,矢量。

球床內(nèi)氣體的動量方程為:

(4)

式中:ui為氣體速度在不同方向上的分量;μ為擴(kuò)散系數(shù);P為流體區(qū)域的壓力;Rxi為對應(yīng)方向上由經(jīng)驗(yàn)關(guān)系式確定的流動阻力;γ為用于確定不同方向是否存在重力加速度的變量。左側(cè)2項(xiàng)分別為非穩(wěn)態(tài)項(xiàng)和對流項(xiàng);右側(cè)5項(xiàng)分別為動量擴(kuò)散項(xiàng)、壓力梯度項(xiàng)、粘性項(xiàng)、由于流固摩擦引起的壓力損失項(xiàng)以及重力項(xiàng)[7]。

結(jié)合準(zhǔn)穩(wěn)態(tài)假設(shè)以及球床壓降公式,球床內(nèi)氣體的動量守恒方程為:

(5)

式中:R為由經(jīng)驗(yàn)關(guān)系式確定的流動阻力,矢量,通過引入球床阻力系數(shù)的經(jīng)驗(yàn)關(guān)系式,可得:

(6)

式中:ψ為球床的摩擦阻力系數(shù);dk為燃料球直徑;ε為球床孔隙率[8]。

氣體溫度場的求解在空間上同樣采用中心差分格式,計(jì)算網(wǎng)格劃分與固體計(jì)算網(wǎng)格劃分一致,待求解的溫度位于材料網(wǎng)格端點(diǎn)。在求解氣體壓力場時,差分網(wǎng)格與材料網(wǎng)格重合,待求解的壓力值位于材料網(wǎng)格的中心處,待求解的質(zhì)量流量位于網(wǎng)格邊界上。

圖2為程序的計(jì)算流程圖,DAYU3D程序包含初始化模塊、固體溫度計(jì)算模塊和對流換熱計(jì)算模塊等,這些模塊的功能相對獨(dú)立,通過程序內(nèi)部的物理量傳遞實(shí)現(xiàn)耦合。在求解固體溫度場時,程序?qū)⒘黧w的溫度場、對流換熱系數(shù)以及功率分布等作為已知量,計(jì)算固體溫度場。在進(jìn)行對流計(jì)算時,同樣將固體溫度場作為已知條件來進(jìn)行計(jì)算。

圖2 DAYU3D計(jì)算流程Fig.2 Flowchart of the DAYU3D calculation process

2 SANA實(shí)驗(yàn)建模

SANA實(shí)驗(yàn)裝置建成于20世紀(jì)90年代,是由德國于利希研究中心建立并用于研究球床高溫氣冷堆內(nèi)的熱量傳遞以及事故工況下的熱傳輸和余熱載出問題[9]。

2.1 SANA實(shí)驗(yàn)裝置介紹

實(shí)驗(yàn)裝置由直徑1.5 m、高1 m的石墨球床、4個加熱元件(1個中心熱元件和3個徑向加熱元件)、頂部和底部絕熱層以及鋼殼構(gòu)成。球床區(qū)域由約9 500個直徑為6 cm的石墨球堆積構(gòu)成,孔隙率約為0.41。3個徑向加熱元件在周向上均勻布置在半徑50 cm處。球床內(nèi)的熱源通過4根加熱元件進(jìn)行提供,中心加熱元件為內(nèi)徑22 mm,外徑32 mm的石墨加熱管,徑向加熱元件為內(nèi)徑10 mm、外徑20 mm的石墨加熱管,加熱元件與球床之間由石墨保護(hù)管隔開。球床布置在一個圓柱形的鋼制容器中,頂部與底部的絕緣材料用以確保大部分熱量水平流經(jīng)球床,在球床內(nèi)部填充有惰性氣體(氦氣或氮?dú)?來模擬球床高溫堆的實(shí)際情況并且避免由高溫引起的石墨腐蝕。實(shí)驗(yàn)裝置如圖3所示,SANA實(shí)驗(yàn)裝置內(nèi)的溫度測量采用熱電偶,分別布置在容器壁、頂部與底部絕緣層、加熱元件保護(hù)管以及球床內(nèi)部的不同高度[10]。

圖3 SANA實(shí)驗(yàn)裝置示意Fig.3 Schematic diagram of the SANA facility

利用該裝置,于利希研究中心開展了在對稱及非對稱工況下的球床加熱實(shí)驗(yàn)。在對稱實(shí)驗(yàn)中,僅通過中心加熱元件對球床進(jìn)行加熱,該工況下的熱源分布滿足二維軸對稱假設(shè),裝置內(nèi)溫度場在周向上保持一致。采用氣體He/N2填充,對稱加熱的實(shí)驗(yàn)工況如表1所示[11]。

表1 SANA對稱加熱實(shí)驗(yàn)Table 1 Symmetric heating experiments by SANA

實(shí)驗(yàn)裝置的結(jié)構(gòu)及溫度測點(diǎn)布置如圖4所示,在軸向高度為9 cm、50 cm以及90 cm分別布置了底部、中部以及頂部溫度測點(diǎn)。

圖4 對稱加熱實(shí)驗(yàn)溫度測點(diǎn)布置Fig.4 Measurement points in the symmetric heating experiment

非對稱加熱實(shí)驗(yàn)中,球床由直徑為60 mm石墨球堆積而成,采用氣體He填充,中心與徑向加熱元件同時對球床進(jìn)行加熱,使得熱源在周向上出現(xiàn)周期性變化,球床內(nèi)的溫度分布也因此呈現(xiàn)出三維效應(yīng),實(shí)驗(yàn)的詳細(xì)工況如表2所示[12]。

表2 SANA非對稱加熱實(shí)驗(yàn)Table 2 Asymmetric heating experiments by SANA

對應(yīng)的溫度測點(diǎn)布置如圖5所示,測點(diǎn)分別位于含徑向加熱元件以及無徑向加熱元件截面軸向高度50 cm處。

圖5 非對稱加熱實(shí)驗(yàn)溫度測點(diǎn)布置Fig.5 Measurement points in the asymmetric heating experiment

2.2 SANA裝置建模

根據(jù)SANA實(shí)驗(yàn)裝置的結(jié)構(gòu)以及幾何參數(shù),建立了如圖6的SANA實(shí)驗(yàn)計(jì)算模型,模型中的各部件幾何尺寸、材料如表3所示。由于實(shí)驗(yàn)所使用的環(huán)形加熱元件的內(nèi)部材料并未在實(shí)驗(yàn)報(bào)告中給出,為了保證計(jì)算建模與實(shí)驗(yàn)裝置結(jié)構(gòu)的一致性,此處將中心加熱元件簡化為半徑16 mm的圓柱形加熱棒。圖6中材料1為孔隙率約0.39的隨機(jī)填充球床,在靠近內(nèi)側(cè)石墨保護(hù)管以及外側(cè)鋼殼處的填充率根據(jù)SANA實(shí)驗(yàn)報(bào)告分別定為0.22和0.52,徑向?qū)挾葹?0 mm,球床的等效導(dǎo)熱系數(shù)由Zehner-Schlünder關(guān)系確定[13]。

表3 SANA裝置計(jì)算模型幾何參數(shù)

注:1.中心球床,2.中心加熱元件,3.保護(hù)管,4.頂絕熱層,5～9.底絕熱層,10.對流邊界,11.內(nèi)側(cè)球床,12.外側(cè)球床,13.電極,14.鋼殼,15.石墨。圖6 對稱加熱實(shí)驗(yàn)計(jì)算模型Fig.6 Model of the symmetric heating experiment

加熱元件與保護(hù)管之間的熱量傳遞主要通過熱輻射進(jìn)行,在DAYU3D程序中,熱輻射模型被納入到了等效導(dǎo)熱系數(shù)的計(jì)算模型中,圖6中材料15為發(fā)射率0.90的石墨材料。實(shí)驗(yàn)裝置頂部及底部絕熱材料的導(dǎo)熱系數(shù)見文獻(xiàn)[11]。

鋼殼與外邊界處的對流換熱系數(shù)根據(jù)SANA實(shí)驗(yàn)報(bào)告[11]以及 Baggemann的研究[13]中給出的推薦值定為18.4 W/(m2·K),外邊界溫度為定溫20 ℃。裝置中加熱元件的實(shí)際熱功率約為額定功率的90%,計(jì)算中采用的熱源密度為加熱元件的實(shí)際熱功率[14]。

對于非對稱加熱工況的計(jì)算模型,無徑向加熱元件柵元的材料網(wǎng)格劃分與圖6一致,包含徑向加熱元件的扇區(qū)的材料網(wǎng)格劃分如圖7所示,在距離裝置中心50 cm處設(shè)置有徑向的加熱元件,程序在三維圓柱坐標(biāo)系下進(jìn)行建模,徑向上的環(huán)形加熱元件在模型中被等效為3個內(nèi)徑48.5 cm,外徑51.5 cm,高度100 cm,θ方向柵元大小為4°的四棱柱加熱區(qū)域,等效過程中保證徑向加熱區(qū)域的功率密度以及換熱面積與徑向加熱元件近似一致。計(jì)算模型在徑向劃分了18個柵元,軸向劃分了35個柵元,θ方向劃分了24個柵元。

注:1.中心球床,2.中心加熱元件,3.保護(hù)管,4.頂絕熱層,5～9.底絕熱層,10.對流邊界,11.內(nèi)側(cè)球床,12.外側(cè)球床,13.電極,14.鋼殼,15.石墨,16.保護(hù)管,17.徑向加熱元件。圖7 非對稱加熱實(shí)驗(yàn)計(jì)算模型Fig.7 Model for the asymmetric heating experiment

3 SANA實(shí)驗(yàn)結(jié)果模擬

3.1 SANA實(shí)驗(yàn)對稱工況模擬

本文首先選擇了大功率水平下的對稱加熱實(shí)驗(yàn)進(jìn)行驗(yàn)證。在大功率水平下,自然對流對于計(jì)算模型的影響可以忽略,同時,更大的溫度梯度有助于對模型存在的問題進(jìn)行修正。因此,此處選擇了30 kW中心加熱、60 mm石墨球床、He作為保護(hù)氣體的對稱工況進(jìn)行了模擬。

圖8為30 kW功率水平下,不同高度處程序計(jì)算結(jié)果與實(shí)驗(yàn)測量值的對比。可以發(fā)現(xiàn),兩者總體上符合較好,計(jì)算誤差在10%以內(nèi),計(jì)算得到的球床內(nèi)部溫度變化趨勢與實(shí)驗(yàn)結(jié)果保持一致。隨著半徑的增大,球床的換熱面積增大,導(dǎo)熱系數(shù)隨溫度的降低而降低,在靠近對流邊界出,球床溫度變化趨于平緩,證明了程序?qū)τ谇虼矁?nèi)部的熱量傳遞的模擬是可靠的。

圖8 SANA裝置30 kW對稱加熱實(shí)驗(yàn)計(jì)算結(jié)果Fig.8 Results for the symmetric heating experiment of the SANA facility at 30 kW

3.2 SANA實(shí)驗(yàn)非對稱工況模擬

為了排除自然對流對于計(jì)算結(jié)果的影響,對于非對稱工況實(shí)驗(yàn),選取了表2中的第4組實(shí)驗(yàn)進(jìn)行了計(jì)算,球床填充的惰性氣體為氦氣。同時,為了與DAYU2D進(jìn)行對比,還將徑向的熱源在θ方向上進(jìn)行平均,得到了均勻化處理之后的二維計(jì)算結(jié)果。圖9為球床中部(Z=50 cm)不同截面處溫度的計(jì)算值與實(shí)驗(yàn)結(jié)果的對比,其中,0°截面和60°截面分別對應(yīng)圖5含徑向加熱元件截面以及無徑向加熱元件截面。

在球床內(nèi)部靠近加熱元件區(qū)域,實(shí)驗(yàn)值略高于計(jì)算結(jié)果,這種差異主要是由于計(jì)算模型將加熱元件近似為圓柱形加熱棒,導(dǎo)致加熱元件的功率密度略低于實(shí)際的功率水平。

從圖9(a)可以發(fā)現(xiàn),0°截面處固體溫度隨著向徑向加熱元件(R=50 cm)靠近,溫度曲線先下降,隨后開始上升,并在經(jīng)過了R=50 cm處后,溫度沿徑向繼續(xù)下降。60°截面處由于遠(yuǎn)離徑向加熱元件,溫度并未出現(xiàn)明顯上升,變化趨勢與對稱工況下的固體溫度分布曲線類似。徑向熱源均勻化處理后計(jì)算得到的溫度曲線介于0°與60°截面的溫度分布曲線之間,與二者的溫度相比均有明顯的差異。

DAYU3D和DAYU2D程序計(jì)算得到的徑向加熱元件處的溫度如表4所示,均勻化處理后計(jì)算得到的最高溫度與實(shí)際溫度823 ℃相差超過150 ℃,表明了三維效應(yīng)對于球床溫度的影響在功率水平較高的工況下不可忽略,二維均勻化處理對三維工況下球床的實(shí)際溫度分布反映不夠充分,處理球床高溫堆中實(shí)際存在的三維熱工問題存在困難。

表4 徑向加熱元件處溫度Table 4 Temperature at the radial heating element

程序?qū)ANA實(shí)驗(yàn)非對稱工況的計(jì)算結(jié)果與測量值符合的較好,溫度曲線的變化趨勢接近,實(shí)驗(yàn)測量值均位于10%誤差線內(nèi)。在0°截面處,加熱棒(R=50 cm)附近的溫度值略低于測量值,這主要是由于模型將徑向加熱元件近似為四棱柱導(dǎo)致計(jì)算得到的功率密度小于加熱元件實(shí)際的功率密度。在靠近側(cè)面對流邊界處,計(jì)算結(jié)果與實(shí)驗(yàn)值的差異主要是由邊界處的對流換熱系數(shù)的不確定性引起的。

4 結(jié)論

1)本文利用國際公開的SANA實(shí)驗(yàn)數(shù)據(jù),驗(yàn)證了DAYU3D程序的三維熱工特性分析能力。球床內(nèi)的非均勻效應(yīng)對內(nèi)部溫度場分布有重要影響。

2)相較于二維熱工分析程序的均勻化處理方式,三維熱工分析程序能夠更精確地描述球床的熱工特性,能更準(zhǔn)確、合理的分析和評價(jià)真實(shí)球床高溫氣冷堆中的三維熱工現(xiàn)象。

在實(shí)現(xiàn)了三維熱工分析程序計(jì)算功能的基礎(chǔ)上,后續(xù)還需開展更深入的工作,包括完善SANA計(jì)算模型加熱元件的近似處理,對程序的功能開展更大范圍的三維算例驗(yàn)證,與三維中子動力學(xué)計(jì)算程序耦合以探究三維熱工程序在球床式高溫氣冷堆熱工設(shè)計(jì)及安全分析中的應(yīng)用等,為進(jìn)一步提高高溫氣冷堆經(jīng)濟(jì)性和安全性提供支持。