基于SANA實驗的三維熱工分析程序計算驗證

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

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

球床式高溫氣冷堆是一種具有第4代核能系統特征的先進反應堆,其固有安全性能夠確保在各種事故工況下,堆芯的燃料最高溫度不會超過設計限值,從而避免大量放射性泄漏[1]。二維系統分析程序THERMIX和TINTE等被廣泛應用于高溫氣冷堆內的熱工分析,但這些程序的建模和計算都是基于二維的圓柱對稱假設,難以處理高溫氣冷堆內存在的由于功率分布不均勻或幾何結構不對稱引起的三維熱工問題。因此,開展針對球床式高溫氣冷堆的三維熱工水力分析十分必要[2]。商用流體動力學計算軟件如FLUENT和CFX等可以應用于球床高溫氣冷堆的局部三維熱工水力分析,但全堆建模的時效性、復雜性以及難以與中子動力學計算程序進行耦合等問題仍然存在。因此,開發能夠實現全堆快速三維熱工計算的系統分析程序對于下一代高溫氣冷堆的設計工作具有重要價值[3]。

二維熱工分析程序DAYU2D在吸收THERMIX軟件包的基礎上開發。基于該程序框架,本文根據堆芯內部的守恒方程和適用于球床式高溫氣冷堆的相關經驗關系式,建立了三維熱工水力計算模型,開發了三維熱工分析程序DAYU3D,并針對德國SANA基準實驗的部分工況進行了計算分析,驗證了程序的計算功能。

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

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

1.1 三維熱傳導計算模型

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

(1)

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

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

圖1 固體節點(I,J,K)處的差分網格Fig.1 Differential grid at solid node (I, J, K)

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

1.2 三維對流計算模型

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

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

(2)

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

氣體的質量守恒方程可以表示為:

(3)

式中G為氣體的質量流密度,矢量。

球床內氣體的動量方程為:

(4)

式中:ui為氣體速度在不同方向上的分量;μ為擴散系數;P為流體區域的壓力;Rxi為對應方向上由經驗關系式確定的流動阻力;γ為用于確定不同方向是否存在重力加速度的變量。左側2項分別為非穩態項和對流項;右側5項分別為動量擴散項、壓力梯度項、粘性項、由于流固摩擦引起的壓力損失項以及重力項[7]。

結合準穩態假設以及球床壓降公式,球床內氣體的動量守恒方程為:

(5)

式中:R為由經驗關系式確定的流動阻力,矢量,通過引入球床阻力系數的經驗關系式,可得:

(6)

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

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

圖2為程序的計算流程圖,DAYU3D程序包含初始化模塊、固體溫度計算模塊和對流換熱計算模塊等,這些模塊的功能相對獨立,通過程序內部的物理量傳遞實現耦合。在求解固體溫度場時,程序將流體的溫度場、對流換熱系數以及功率分布等作為已知量,計算固體溫度場。在進行對流計算時,同樣將固體溫度場作為已知條件來進行計算。

圖2 DAYU3D計算流程Fig.2 Flowchart of the DAYU3D calculation process

2 SANA實驗建模

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

2.1 SANA實驗裝置介紹

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

圖3 SANA實驗裝置示意Fig.3 Schematic diagram of the SANA facility

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

表1 SANA對稱加熱實驗Table 1 Symmetric heating experiments by SANA

實驗裝置的結構及溫度測點布置如圖4所示,在軸向高度為9 cm、50 cm以及90 cm分別布置了底部、中部以及頂部溫度測點。

圖4 對稱加熱實驗溫度測點布置Fig.4 Measurement points in the symmetric heating experiment

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

表2 SANA非對稱加熱實驗Table 2 Asymmetric heating experiments by SANA

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

圖5 非對稱加熱實驗溫度測點布置Fig.5 Measurement points in the asymmetric heating experiment

2.2 SANA裝置建模

根據SANA實驗裝置的結構以及幾何參數,建立了如圖6的SANA實驗計算模型,模型中的各部件幾何尺寸、材料如表3所示。由于實驗所使用的環形加熱元件的內部材料并未在實驗報告中給出,為了保證計算建模與實驗裝置結構的一致性,此處將中心加熱元件簡化為半徑16 mm的圓柱形加熱棒。圖6中材料1為孔隙率約0.39的隨機填充球床,在靠近內側石墨保護管以及外側鋼殼處的填充率根據SANA實驗報告分別定為0.22和0.52,徑向寬度為30 mm,球床的等效導熱系數由Zehner-Schlünder關系確定[13]。

表3 SANA裝置計算模型幾何參數

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

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

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

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

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

3 SANA實驗結果模擬

3.1 SANA實驗對稱工況模擬

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

圖8為30 kW功率水平下,不同高度處程序計算結果與實驗測量值的對比。可以發現,兩者總體上符合較好,計算誤差在10%以內,計算得到的球床內部溫度變化趨勢與實驗結果保持一致。隨著半徑的增大,球床的換熱面積增大,導熱系數隨溫度的降低而降低,在靠近對流邊界出,球床溫度變化趨于平緩,證明了程序對于球床內部的熱量傳遞的模擬是可靠的。

圖8 SANA裝置30 kW對稱加熱實驗計算結果Fig.8 Results for the symmetric heating experiment of the SANA facility at 30 kW

3.2 SANA實驗非對稱工況模擬

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

在球床內部靠近加熱元件區域,實驗值略高于計算結果,這種差異主要是由于計算模型將加熱元件近似為圓柱形加熱棒,導致加熱元件的功率密度略低于實際的功率水平。

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

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

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

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

4 結論

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

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

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