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

(哈爾濱工程大學(xué) 核安全與仿真技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，黑龍江 哈爾濱 150001)

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

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

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

1 計(jì)算幾何模型

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

圖2 普羅米修斯計(jì)劃空間反應(yīng)堆結(jié)構(gòu)

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

圖3 燃料棒結(jié)構(gòu)

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

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

2 計(jì)算模型及條件設(shè)置

2.1 邊界條件設(shè)置

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

圖4 反應(yīng)堆堆芯功率分布

2.2 網(wǎng)格模型

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

2.3 物理模型

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

圖5 網(wǎng)格結(jié)構(gòu)

3 計(jì)算結(jié)果及分析

3.1 網(wǎng)格無(wú)關(guān)性驗(yàn)證

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

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

3.2 計(jì)算模型驗(yàn)證

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

表2 計(jì)算模型驗(yàn)證算例參數(shù)

圖7 仿真驗(yàn)證算例與文獻(xiàn)實(shí)驗(yàn)數(shù)據(jù)對(duì)比

3.3 堆芯流場(chǎng)分析

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

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

圖9 出口管橫截面速度場(chǎng)

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

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

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

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

圖11 堆芯溫度場(chǎng)分布

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

圖12 堆芯壓強(qiáng)分布

堆芯入口處，由于入口管道幾乎垂直于內(nèi)壓力容器，因此入流冷卻劑直接撞擊內(nèi)壓力容器，產(chǎn)生較大能量損失。從優(yōu)化的角度考慮，可改變?nèi)肟诠艿婪较颍蛊湎蛏戏絻A斜，與壓力容器下降環(huán)段呈鈍角，以避免入流冷卻劑直接撞擊內(nèi)壓力容器，進(jìn)而減少不必要的能量損失。

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

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

圖13 不同冷卻劑通道內(nèi)溫度分布

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

圖14 冷卻劑通道示意圖

3.4 不確定性分析

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

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

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

4 結(jié)論

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

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

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

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