孟 濤，趙富龍，程 坤，曾 陳，譚思超

(哈爾濱工程大學 核安全與仿真技術國防重點學科實驗室，黑龍江 哈爾濱 150001)

近年來，諸如空間反應堆等特種反應堆開始受到越來越多的關注。此類反應堆因功率較小、結構相對簡單，利于反應堆的小型化，因而在對空間體積和重量要求較為嚴格的空間環(huán)境、船舶環(huán)境中具有更為廣闊的應用前景。其中直接布雷頓循環(huán)氣冷反應堆因結構緊湊、體積重量較小、功率適用范圍大等特點而多應用于大功率空間系統(tǒng)[1]。

美俄均提出過多個空間氣冷反應堆研究計劃[2-4]。2003年3月，美國航空航天局簽署了普羅米修斯核能系統(tǒng)研究計劃[1]，該計劃作為木星冰衛(wèi)星探測器計劃的后續(xù)計劃，目標在于研制一款安全、可靠、能執(zhí)行多種外太陽系探索任務的深空探測器。該探測器以核電推進為動力，可攜帶任務載荷不小于1.5 t，并要求任務相關技術可應用到月球與火星表面供電任務以及空間貨物運輸任務上。雖然該任務開展僅兩年后便被終止，但該任務對空間反應堆相關技術進行了廣泛的調研，針對空間用反應堆堆芯類型、熱電轉換方式、反應堆結構材料等，給出了全面、詳細的技術調研報告，能為相關研究提供很好的幫助。另一方面俄羅斯也提出了兆瓦級空間核反應堆計劃[5]，該計劃旨在開發(fā)兆瓦級功率空間供電與推進反應堆系統(tǒng)，能為空間探索及輸送任務提供動力。

為揭示空間氣冷反應堆堆芯流動換熱機理，本文基于美國普羅米修斯計劃的公開信息[1,5]，建立相應的空間反應堆堆芯CFD模型，針對空間反應堆堆芯流動與傳熱特性開展計算分析，評估現(xiàn)有設計中仍待優(yōu)化之處，并提出相關的優(yōu)化建議，為未來空間氣冷反應堆設計提供參考。

1 計算幾何模型

普羅米修斯計劃提出了多個氣冷堆方案，包括基準方案、低功率方案、開放流通區(qū)域方案等，但計劃終止時并未給出最終設計方案[6]。為分析空間反應堆的流動換熱特性，本文參照基準設計方案進行數(shù)值仿真計算。反應堆堆芯內燃料棒以正六邊形排布，由內至外共排布有7圈燃料棒，最外層1圈燃料棒在正六邊形的6個角各缺少1根燃料棒，燃料排布如圖1所示。中心區(qū)域不排布燃料棒，僅設置1根控制棒，且控制棒通道與冷卻劑通道隔離，堆芯幾何模型如圖2所示。

圖2 普羅米修斯計劃空間反應堆結構

反應堆由內、外兩層壓力容器構成，內外兩層壓力容器通過開有孔洞的柵板連接。燃料棒結構及冷卻劑流動方式如圖3所示。每根燃料棒具有獨立的環(huán)形冷卻劑通道，燃料棒通過入口、出口兩個端部與環(huán)形冷卻劑通道及堆芯上下支撐板固定。堆芯主要參數(shù)信息列于表1。

圖3 燃料棒結構

參數(shù)數(shù)值燃料棒外徑,mm19.65冷卻劑通道厚度,mm2.16燃料棒長度,m1.118燃料棒芯塊直徑,mm18.19燃料活性區(qū)長度,m0.608燃料棒數(shù)量288堆芯壓力容器外徑,mm618.1壓力容器內環(huán)隙厚度,m0.3冷卻劑質量流量,kg/s6.16堆芯入口溫度,K911堆芯出口溫度,K1150反應堆熱功率,kW1002

由于文獻中的堆芯幾何參數(shù)信息有限，因此本文使用的部分堆芯結構參數(shù)為推測數(shù)據(jù)，包括內壓力容器內徑、冷卻劑主管道內徑和燃料棒出口端結構尺寸。數(shù)據(jù)的推測依據(jù)為普羅米修斯計劃反應堆結構的幾何比例。推測數(shù)據(jù)對堆芯內部流動與換熱影響較小，主要影響進出口處的流速及流場特征。

2 計算模型及條件設置

2.1 邊界條件設置

由于該反應堆的冷卻劑通道為獨立通道，故選取單根燃料棒及其冷卻劑通道作為單通道模型進行計算。首先建立了單通道模型進行仿真，以進行網(wǎng)格無關性驗證。入口條件均設置為質量流量入口，經(jīng)過計算確定單根燃料棒通道平均冷卻劑流量為0.021 389 kg/s，相應的1/6堆芯冷卻劑流量為1.026 67 kg/s，入口溫度為911 K，出口條件設置為壓力出口，壓強為1 898 kPa。單通道燃料棒加熱方式為余弦加熱，積分加熱功率為3.479 kW。1/6堆芯的加熱方式采用真實堆芯功率分布，堆芯出口平均溫度控制在1 150 K。堆芯功率分布通過使用SuperMC蒙特卡羅軟件[7-8]計算得到，如圖4所示。在熱工水力計算中使用線性表格插值的方式調用功率分布。計算中使用He-Xe混合氣體作為冷卻劑[9-11]，混合摩爾質量參照普羅米修斯方案選定為31.5 g/mol。

圖4 反應堆堆芯功率分布

2.2 網(wǎng)格模型

多面體網(wǎng)格與其他非結構網(wǎng)格相比具有最少的網(wǎng)格數(shù)量、較少的計算時間和較高的計算精度，因此本文采用多面體網(wǎng)格，網(wǎng)格劃分采用自適應加密方法，可根據(jù)結構特點自動加密網(wǎng)格，基準網(wǎng)格尺寸為6 mm，最小網(wǎng)格尺寸為基準網(wǎng)格尺寸的25%，近壁面區(qū)采用邊界層網(wǎng)格，邊界層網(wǎng)格層數(shù)為5，總厚度為1.2 mm，最終劃分的局部網(wǎng)格結構如圖5所示。

2.3 物理模型

由于目前關于He-Xe混合氣體的相關實驗公開文獻較少[12-16]，且實驗精度難以保證，因此，本文使用數(shù)值計算方法模擬空間氣冷反應堆堆芯冷卻劑流動情況，并參照文獻[15]中的對流換熱實驗數(shù)據(jù)，進而獲得堆芯流動換熱特性，為空間反應堆的優(yōu)化設計提供指導。本文使用適用范圍較廣的可實現(xiàn)k-ε模型進行數(shù)值計算?？蓪崿F(xiàn)k-ε模型是兩方程雷諾平均湍流模型，由Shih等[17]基于標準k-ε模型推導得出。由于該模型對湍流黏性系數(shù)的計算采用了時均應變率、旋度及旋轉坐標系角速度的函數(shù)，因而相比于標準k-ε模型通常具有更好的計算精度和網(wǎng)格適應性。

圖5 網(wǎng)格結構

3 計算結果及分析

3.1 網(wǎng)格無關性驗證

圖6 參數(shù)隨網(wǎng)格數(shù)量的變化

為確保計算結果的精度和可靠性，首先需進行網(wǎng)格無關性驗證。圖6示出了利用單通道模型計算得到的網(wǎng)格敏感性分析結果。結果表明，當網(wǎng)格數(shù)量大于30萬時，冷卻劑出口最高溫度與平均速度的變化幅度開始顯著變小，相對誤差在3%以內。因此，對于單通道，網(wǎng)格數(shù)量應控制在30萬以上。而對于1/6堆芯，由于其包含54根燃料通道、上下腔室以及下降段，因此本文選擇的網(wǎng)格數(shù)量為2 558萬，網(wǎng)格質量大于0.9。

3.2 計算模型驗證

為確保計算模型的可靠性，需進行計算模型驗證。由于使用He-Xe混合氣體作為冷卻劑的流動換熱特性及仿真屬于細分研究領域，且實驗成本高昂，目前為止公開發(fā)表的實驗文獻及仿真研究極其有限[12-16]。本文參照實驗參數(shù)相對全面的文獻[15]進行計算模型驗證。計算使用可實現(xiàn)k-ε模型，網(wǎng)格采用多面體網(wǎng)格，網(wǎng)格數(shù)量為52.5萬，計算邊界條件列于表2，計算結果示于圖7。由計算結果發(fā)現(xiàn)，冷卻劑中心溫度與文獻[15]中實驗數(shù)據(jù)符合較好，最大溫度誤差小于5%。但加熱管壁面溫度相差較多。通過在流固交界面增加接觸熱阻0.005 m2·K/W后，加熱管壁面溫度在整體趨勢上與文獻[15]中的實驗數(shù)據(jù)能較好地符合。由此說明，使用可實現(xiàn)k-ε模型可較好地模擬He-Xe混合氣體管內流動的流體特性，但仿真得到的壁面溫度會偏低，進而導致計算得到的努塞爾數(shù)Nu或對流換熱系數(shù)h偏低。增加接觸熱阻后，可對壁面溫度進行較好的仿真模擬。由于本文進行的1/6堆芯計算，未考慮流固耦合，故未在計算中增加接觸熱阻。

表2 計算模型驗證算例參數(shù)

圖7 仿真驗證算例與文獻實驗數(shù)據(jù)對比

3.3 堆芯流場分析

為簡化計算，忽略燃料棒、控制棒、格架等結構部件，僅針對流體區(qū)域進行數(shù)值計算。通過計算得到的冷卻劑在入口段和出口段的速度云圖如圖8所示。由圖8a可知，在入口段，流體進入壓力容器間環(huán)隙后直接撞擊在內層壓力容器壁面，形成分流，并在撞擊點附近形成低流速區(qū)。向上方弧形段流動的冷卻劑流速快速降低，向下方流動的冷卻劑通過開孔柵板時由于流通面積的減少而出現(xiàn)加速流動。通過開孔柵板后流速逐漸降低，流動到下方開孔柵板時又出現(xiàn)加速流動，之后流速再逐漸降低，并開始向堆芯中心側偏轉，如圖8b所示。在偏轉流動過程中，靠近弧形壁面的一部分入流冷卻劑繼續(xù)沿壁面向下流動，直至流動到堆芯縱向中軸線位置處，在受到中心控制棒通道壁面阻擋后改變流動方向，沿控制棒通道壁面豎直向上運動，與入流冷卻劑混合后共同流入燃料棒冷卻劑通道。流出燃料棒通道后，不同冷卻劑通道的流體匯聚后通過出口管流出反應堆。

圖8 冷卻劑入口段(a)和出口段(b)速度場

圖9 出口管橫截面速度場

為了解冷卻劑在出口位置處的速度分布情況，圖9示出了出口管橫截面的速度云圖。由圖9可知，在冷卻劑流出管橫截面上，速度分布不均勻，上方速度較高，下方速度較低。

圖10示出了不同高度位置處的堆芯橫截面速度分布結果。在入口管中間位置(圖10a)處，堆芯速度場沿出口管中間軸線呈對稱分布。在圓周方向中間位置處，冷卻劑流速較高，而靠近內壓力容器壁面處的外圍冷卻劑流速較低。冷卻劑流出燃料棒區(qū)域后快速進入出口管，并未在上封頭內大范圍流動，因而使得流速分布不均勻，產(chǎn)生了局部高流速區(qū)和低流速區(qū)。在入口段下方約8 cm位置(圖10b)處，環(huán)狀下降段內冷卻劑流速分布并不均勻，入口管內的冷卻劑進入后受壓力容器壁面的影響形成分流，導致在徑向方向上冷卻劑速度分布不均勻，進而影響沿環(huán)狀下降段流動的冷卻劑的速度分布。

為了解堆芯通道中的溫度分布情況，圖11示出了冷卻劑的溫度分布云圖。由圖11可知，當冷卻劑進入堆芯后，由于徑向加熱的不均勻性，堆芯徑向中間位置處冷卻劑具有最高溫度，而半徑較小和較大位置處的冷卻劑溫度都相對較低。當冷卻劑流入出口端下封頭時，冷卻劑并未發(fā)生明顯攪混，進而造成該位置處冷卻劑溫度分布不均勻。

圖10 入口管段中間位置(a)及下方8 cm處(b)冷卻劑速度場分布

圖11 堆芯溫度場分布

此外，為了解冷卻劑在堆芯通道中流動過程中流動阻力的變化情況，為結構優(yōu)化設計提供參考，圖12示出了堆芯通道內的壓強分布云圖。由圖12可知，冷卻劑在流動過程中的壓降主要發(fā)生在3個位置，1個是環(huán)狀下降段的2處支撐處，另1個則是在堆芯燃料區(qū)域的冷卻劑通道中。

圖12 堆芯壓強分布

堆芯入口處，由于入口管道幾乎垂直于內壓力容器，因此入流冷卻劑直接撞擊內壓力容器，產(chǎn)生較大能量損失。從優(yōu)化的角度考慮，可改變入口管道方向，使其向上方傾斜，與壓力容器下降環(huán)段呈鈍角，以避免入流冷卻劑直接撞擊內壓力容器，進而減少不必要的能量損失。

另一方面，由圖8a可知，冷卻劑流入壓力容器后，一部分向上方的弧段環(huán)隙流動。但弧段環(huán)隙內大部分流體流動較為緩慢。因此，可考慮對上部弧形封頭進行結構優(yōu)化設計，減輕系統(tǒng)質量。

此外，在冷卻劑出口段，冷卻劑攪混程度較輕，導致出口段管道、弧形封頭內存在較大的溫差與速度差，不利于平衡封頭結構和出口管道的內部熱應力，既有可能降低結構壽命，同時還會影響氣輪機工作性能。因此，從優(yōu)化的角度考慮，可考慮設置加強攪混結構或優(yōu)化出口管排布位置，來達到降低冷卻劑出口溫度分布不均勻的目的。具體優(yōu)化方式仍有待進一步計算分析。

圖13 不同冷卻劑通道內溫度分布

為評價不同冷卻劑通道內流體的流動特征，圖13示出了不同冷卻劑通道內冷卻劑溫度隨堆芯高度的變化，圖中標號如圖14所示。由圖13可知，不同通道的溫度變化均呈現(xiàn)相同的升溫趨勢，升溫速率均是先小后大，然后再逐漸減小。越靠近半徑中間位置，加熱功率越大，通道出口溫度越高。計算結果顯示，不同通道出口溫度最大溫差可達76 K。因此，當堆芯出口處攪混不夠充分時，堆芯出口管內冷卻劑將產(chǎn)生較大的溫度不均勻分布，進而影響氣輪機設備的正常運行。

圖14 冷卻劑通道示意圖

3.4 不確定性分析

本文參照文獻[18]中的不確定分析方法，針對3種不同網(wǎng)格情況進行仿真結果的不確定性分析，網(wǎng)格數(shù)量分別為2 558萬(網(wǎng)格1)、4 284萬(網(wǎng)格2)和4 733萬(網(wǎng)格3)。計算得到的冷卻劑出口溫度、平均流速和壓降列于表3。根據(jù)計算結果得到的網(wǎng)格收斂系數(shù)(GCI)列于表4，其中網(wǎng)格1/2表示利用網(wǎng)格1與網(wǎng)格2計算結果得到網(wǎng)格收斂系數(shù)，依次類推。將最大網(wǎng)格收斂系數(shù)除以2即可得到數(shù)值不確定度。由此可知，仿真計算的冷卻劑出口溫度不確定度為0.908 5 K，冷卻劑出口流速不確定度為0.182 m/s，反應堆壓降不確定度為1 419 Pa。

表3 不同網(wǎng)格數(shù)量下的計算結果

表4 網(wǎng)格收斂系數(shù)

4 結論

本文針對美國普羅米修斯計劃提出的反應堆方案建立了三維模型，并使用可實現(xiàn)k-ε模型進行了1/6堆芯的流動換熱計算，堆芯加熱功率使用堆芯穩(wěn)態(tài)物理計算得到的功率分布曲線，通過計算得到了堆芯溫度場、速度場和壓力場的分布情況，并分析了空間氣冷反應堆堆芯的流動換熱特性，得到以下主要結論。

1) 在入口段，流體直接沖擊內壓力容器壁面，形成分流，并在撞擊點附近形成低流速區(qū)。從優(yōu)化的角度考慮，可改變入口管道方向，使其向上方傾斜，與壓力容器下降環(huán)段呈鈍角，以避免入流冷卻劑直接撞擊內壓力容器，進而減少不必要的能量損失。

2) 在堆芯下降段及燃料區(qū)域，冷卻劑速度場沿出口管中間軸線呈對稱分布。上腔室內圓周方向中間位置處冷卻劑流速較高，而兩側靠近內壓力容器壁面處的冷卻劑流速較低，環(huán)狀下降段內冷卻劑流速分布并不均勻，入口管內的冷卻劑進入后受壓力容器壁面的影響形成分流，導致在徑向方向上冷卻劑速度分布不均勻，進而影響沿環(huán)狀下降段流動的冷卻劑的速度分布。

3) 堆芯出口管截面上溫度及速度分布不均勻，溫差可達76 K，從優(yōu)化的角度考慮，可設置加強攪混結構或優(yōu)化出口管排布位置，降低冷卻劑對氣輪機等設備的影響。