基于KRUSTY帶核實(shí)驗(yàn)的熱管冷卻反應(yīng)堆瞬態(tài)分析程序驗(yàn)證

吳攀 歐陽(yáng)澤宇 朱煜 單建強(qiáng)

（西安交通大學(xué) 核安全與運(yùn)行實(shí)驗(yàn)室 西安 710049）

核工業(yè)界目前正在研發(fā)多種微型反應(yīng)堆方案（簡(jiǎn)稱“微堆”），微堆是指系統(tǒng)非常簡(jiǎn)單緊湊，能產(chǎn)生1~20 MW的熱能直接供熱或轉(zhuǎn)換為電能利用的反應(yīng)堆系統(tǒng)，具有模塊化、可運(yùn)輸、自適應(yīng)調(diào)節(jié)等核心關(guān)鍵特征，相較于傳統(tǒng)反應(yīng)堆具有微小便攜、設(shè)計(jì)簡(jiǎn)單、可快速現(xiàn)場(chǎng)部署安裝等優(yōu)勢(shì)，能夠很好地解決國(guó)家的特種能源需求。在眾多微堆類型中，熱管冷卻反應(yīng)堆作為能量密度高、適應(yīng)能力強(qiáng)的核電源系統(tǒng)，能在空間、深海探索等任務(wù)中發(fā)揮重要作用，因而越來(lái)越受到各國(guó)重視［1］。使用高溫?zé)峁芾鋮s堆芯的微型熱管冷卻反應(yīng)堆具有高溫運(yùn)行、系統(tǒng)緊湊、固有安全性高的優(yōu)點(diǎn)，便于模塊化設(shè)計(jì)，可避免單點(diǎn)失效。此外，由于省卻了一回路泵、閥門及管道系統(tǒng)，簡(jiǎn)化了系統(tǒng)布置并節(jié)省了系統(tǒng)體積，可避免傳統(tǒng)商用反應(yīng)堆可能發(fā)生的失水與失流事故［2］。高溫?zé)峁芡ǔ２捎娩?、鈉、鉀等堿金屬作為工質(zhì)。高溫?zé)峁艿倪\(yùn)行性能對(duì)高壓堆的啟動(dòng)、穩(wěn)態(tài)、瞬態(tài)和停堆都有顯著的影響［3］，相關(guān)研究表明，微型熱管冷卻反應(yīng)堆在技術(shù)上較其他微堆概念更成熟，如美國(guó)洛斯阿莫斯國(guó)家實(shí)驗(yàn)室研發(fā)的MegaPower熱管反應(yīng)堆［4］預(yù)計(jì)于2025年進(jìn)行原型驗(yàn)證，Kilopower熱管空間堆［5］已于2018年完成了地面驗(yàn)證測(cè)試。因此，微型熱管冷卻反應(yīng)堆成為了先進(jìn)小型反應(yīng)堆中最具應(yīng)用前景的技術(shù)。

對(duì)于熱管冷卻反應(yīng)堆，通過(guò)熱管的熱傳遞是堆芯裂變熱輸出的唯一方法［6］。為分析微型熱管冷卻反應(yīng)堆在正常和非正常運(yùn)行條件下的系統(tǒng)瞬態(tài)響應(yīng)與安全特性，開展系統(tǒng)安全評(píng)估，需開發(fā)系統(tǒng)瞬態(tài)分析模型和程序［7］?，F(xiàn)已有針對(duì)熱管冷卻雙模式空間堆（HP-BSNR）開發(fā)的熱工水力分析程序STHA_HPBSNR［8］、針對(duì)熱管分段式熱電模塊轉(zhuǎn)換器（Heat Pipes Segmented Thermoelectric Module Converters，HP-STMCs）空間堆開發(fā)了熱工水力計(jì)算程序［9］、針對(duì)熱管冷卻空間反應(yīng)堆電源系統(tǒng)（Scalable Alkali metal Thermal-To-Electric Conversion Integrated Reactor Space Power Systems，SAIRS）空間堆開發(fā)的系統(tǒng)瞬態(tài)分析程序TAPIRS［10］等，但由于國(guó)內(nèi)外熱管冷卻反應(yīng)堆相關(guān)實(shí)驗(yàn)開展較晚，這些程序僅對(duì)反應(yīng)堆方案的穩(wěn)態(tài)設(shè)計(jì)值進(jìn)行了初步驗(yàn)算，未經(jīng)過(guò)真實(shí)實(shí)驗(yàn)數(shù)據(jù)的驗(yàn)證。千瓦級(jí)熱管反應(yīng)堆實(shí)驗(yàn)（Kilowatt Reactor Using Stirling Technology，KRUSTY）是美國(guó)國(guó)家航空航天局（National Aeronautics and Space Administration，NASA）針對(duì)1~10 kWe級(jí)小型空間核反應(yīng)堆電源Kilopower項(xiàng)目開展的實(shí)驗(yàn)項(xiàng)目，該實(shí)驗(yàn)是當(dāng)前世界上唯一公開發(fā)表的基于斯特林技術(shù)的熱管冷卻反應(yīng)堆帶核實(shí)驗(yàn)，其結(jié)構(gòu)設(shè)計(jì)如圖1所示。目前，KRUSTY熱管堆系統(tǒng)的28 h帶核實(shí)驗(yàn)數(shù)據(jù)已經(jīng)公開發(fā)布［11-12］，為熱管冷卻反應(yīng)堆系統(tǒng)分析程序的驗(yàn)證提供了十分有力的參考。

針對(duì)當(dāng)前熱管堆系統(tǒng)分析程序驗(yàn)證工作不足的問(wèn)題，本文中介紹了自主研制的熱管冷卻反應(yīng)堆系統(tǒng)瞬態(tài)分析程序TAPIRS-D的關(guān)鍵理論模型，并采用TAPIRS-D對(duì)KRUSTY熱管堆系統(tǒng)的堆芯傳熱系統(tǒng)和熱管傳熱系統(tǒng)進(jìn)行了詳細(xì)建模，模擬計(jì)算了帶核實(shí)驗(yàn)中的多個(gè)運(yùn)行工況，包括熱管堆冷啟動(dòng)工況、負(fù)載變化瞬態(tài)工況、熱管失效事故工況、反應(yīng)性引入事故工況和熱阱喪失事故工況，使用文獻(xiàn)［11-12］中的實(shí)驗(yàn)數(shù)據(jù)評(píng)估了TAPIRS-D程序的堆芯功率瞬變模型、堆芯瞬態(tài)傳熱模型以及熱管瞬態(tài)傳熱模型等的正確性。

1 TAPIRS-D程序關(guān)鍵模型介紹

TAPIRS-D程序是基于TAPIRS空間堆系統(tǒng)瞬態(tài)分析程序［10］進(jìn)行模型改進(jìn)優(yōu)化二次開發(fā)的熱管堆系統(tǒng)分析程序，可針對(duì)采用多種熱電轉(zhuǎn)換方式（包括靜態(tài)和動(dòng)態(tài)熱電轉(zhuǎn)化）、堆芯布置方式和廢熱排出方式（水冷或者輻射散熱）的熱管堆系統(tǒng)進(jìn)行建模計(jì)算。TAPIRS-D程序的模型框架如圖2所示，程序的分析模型主要包括堆芯功率瞬變模型、堆芯傳熱模型、熱管傳熱模型、熱電轉(zhuǎn)換模型及廢熱排出模型，對(duì)于不同模塊建立的微分方程組使用四階標(biāo)準(zhǔn)龍格庫(kù)塔方法進(jìn)行聯(lián)合求解。

圖2 TAPIRS-D程序模型框架Fig.2 Model framework of TAPIRS-D program

由于KRUSTY熱管冷卻反應(yīng)堆實(shí)驗(yàn)的模擬主要涉及的模塊包括功率計(jì)算模塊、堆芯傳熱模塊以及熱管傳熱模塊，本節(jié)將簡(jiǎn)要介紹TAPIRS-D程序在這三個(gè)方面的基本理論模型。

反應(yīng)堆的堆芯功率瞬變模型采用點(diǎn)堆中子動(dòng)力學(xué)模型，使用同時(shí)適用于剛性顯著情況與剛性不顯著情況的高次端點(diǎn)浮動(dòng)法［13］用于求解考慮了6組緩發(fā)中子的點(diǎn)堆中子動(dòng)力學(xué)方程，堆芯的反應(yīng)性反饋考慮燃料與熱管等材料的溫度反應(yīng)性反饋。堆芯傳熱模型采用基于傅里葉定律和能量守恒定律的二維溫度控制方程［14］，使用顯式差分方法可得到關(guān)于各控制體節(jié)點(diǎn)的溫度微分方程組。熱管傳熱模型則通過(guò)熱管在凝固態(tài)啟動(dòng)階段的自由分子流模型、“平面前鋒”模型和正常工作階段的熱阻網(wǎng)格模型，全面模擬熱管從啟動(dòng)到正常運(yùn)行的全過(guò)程傳熱性能。

1.1 堆芯功率瞬變模型

1.1.1 點(diǎn)堆中子動(dòng)力學(xué)模型

核反應(yīng)堆的反應(yīng)性快速變化可能會(huì)導(dǎo)致堆芯功率的快速變化，核反應(yīng)堆裂變功率的控制是利用緩發(fā)中子的作用實(shí)現(xiàn)的［15］。由于微型堆的體積小，反應(yīng)堆的中子動(dòng)力學(xué)特性的空間效應(yīng)并不明顯，同時(shí)為了實(shí)現(xiàn)一維的瞬態(tài)系統(tǒng)分析程序的快速計(jì)算，TAPIRS-D采用考慮6組緩發(fā)中子的點(diǎn)堆中子動(dòng)力學(xué)模型來(lái)計(jì)算反應(yīng)堆功率大小及瞬態(tài)變化：

式中：N（t）為中子密度（功率水平），cm-3；ρ（t）為反應(yīng)性，$；β為緩發(fā)中子總份額；λi為6組緩發(fā)中子先驅(qū)核的衰變常數(shù)，s-1；Ci為第i組緩發(fā)中子先驅(qū)核的濃度，cm-3；Λ為瞬發(fā)中子的平均壽命，s；βi為6組緩發(fā)中子份額。

1.1.2 反應(yīng)性反饋模型

反應(yīng)性反饋模型主要考慮堆芯燃料與基體材料的溫度反應(yīng)性反饋。在快堆運(yùn)行期間，堆芯大小會(huì)發(fā)生變化，這主要是因?yàn)槎褍?nèi)溫度變化引起構(gòu)件的膨脹或移動(dòng)，從而引起堆芯各部分中子吸收能力發(fā)生變化，考慮溫度反應(yīng)性反饋的堆芯總反應(yīng)性計(jì)算式如下：

式中：αF為燃料溫度反應(yīng)性反饋系數(shù)，$·K-1為燃料平均溫度，K；為燃料基準(zhǔn)溫度；αC為基體溫度反應(yīng)性反饋系數(shù)，$·K-1為基體平均溫度，K為基體基準(zhǔn)溫度，K；ρext為外部引入的反應(yīng)性，$。

1.1.3 控制鼓調(diào)節(jié)模型

熱管反應(yīng)堆常用的外部反應(yīng)性控制方法是使用步進(jìn)電機(jī)控制的控制鼓系統(tǒng)，控制鼓的主體是轉(zhuǎn)軸可在0o~180o內(nèi)旋轉(zhuǎn)的圓柱形氧化鈹（BeO）反射層轉(zhuǎn)鼓，側(cè)面一定扇面區(qū)域附有碳化硼（B4C）吸收體。當(dāng)控制鼓朝外正向轉(zhuǎn)動(dòng)時(shí)，吸收體遠(yuǎn)離堆芯，堆芯中子更多地被反射，因此會(huì)向堆芯引入正反應(yīng)性，反之則會(huì)引入負(fù)反應(yīng)性。

控制鼓引入的外部反應(yīng)性價(jià)值可以擬合為與控制鼓轉(zhuǎn)角大小相關(guān)的多項(xiàng)式函數(shù)，圖3所示為一典型控制鼓設(shè)計(jì)的反應(yīng)性價(jià)值ρCD與轉(zhuǎn)角θ間關(guān)系曲線圖，其擬合關(guān)系式的形式如下：

圖3 控制鼓反應(yīng)性價(jià)值與轉(zhuǎn)角之間的關(guān)系曲線Fig.3 Relationship curve between the reactivity value of control drum and the turning angle

1.2 堆芯傳熱模型

熱管堆的堆芯結(jié)構(gòu)相對(duì)復(fù)雜，大型熱管堆的堆芯一般由一定比例的燃料與熱管構(gòu)成的燃料組件布置而成，堆芯放置在隔離層中，外層布置有BeO等中子反射材料組成的控制鼓。由于堆芯存在著徑向?qū)幔覍?shí)驗(yàn)中反射層與環(huán)境的熱損失一般不可忽略，在失去廢熱排出系統(tǒng)的情況下，堆芯余熱也是通過(guò)徑向反射層導(dǎo)出堆芯的，因此首先建立了堆芯徑向傳熱模型，將六邊形堆芯區(qū)域等效為圓形，并將整個(gè)堆芯區(qū)域劃分為熱區(qū)、平均區(qū)、冷區(qū)、容器壁、反射層五個(gè)等效計(jì)算區(qū)域（圖4），每個(gè)計(jì)算區(qū)域內(nèi)可包含燃料棒、熱管、基體、BeO棒等組件。計(jì)算區(qū)域的一維徑向?qū)峥刂品匠掏ㄟ^(guò)傅里葉定律和能量守恒定律得到：

圖4 堆芯徑向分區(qū)示意Fig.4 Diagram of radial partitioning of the core

式中：V為計(jì)算區(qū)域體積，m3為計(jì)算區(qū)域平均溫度，K；r為區(qū)域等效半徑，m；S為內(nèi)熱源項(xiàng)，包括燃料產(chǎn)熱、熱管傳熱、對(duì)流散熱功率等，W。

1.3 熱管傳熱模型

反應(yīng)堆堆芯的熱量由熱管導(dǎo)出，熱管是一個(gè)比較復(fù)雜的系統(tǒng)，其工作全過(guò)程包含了許多傳熱模式，比如自由分子流換熱模式、啟動(dòng)模式和正常換熱模式等。高溫堿金屬熱管運(yùn)行在700~1 600 K的溫度范圍內(nèi)，其內(nèi)部壓力小于0.1 MPa。熱管的傳熱計(jì)算需要涵蓋從冷態(tài)啟動(dòng)到正常運(yùn)行的所有過(guò)程，根據(jù)工質(zhì)轉(zhuǎn)變溫度Ttr，熱管傳熱模型包括自凝固態(tài)啟動(dòng)階段的自由分子流模型、“平面前鋒”模型和正常工作階段的熱阻網(wǎng)格模型。對(duì)于高溫?zé)峁軉?dòng)早期的自由分子流階段，管內(nèi)以蒸發(fā)段工熔火質(zhì)的熔化升溫過(guò)程為主，蒸汽稀薄，絕熱段、冷凝段工質(zhì)仍處于凝固狀態(tài)，因此該階段可使用能量守恒方程描述：

式中：C(Te)為熱管溫度等于Te時(shí)的單位長(zhǎng)度熱容量，J·（m·K）-1；Te為蒸發(fā)段平均溫度，K；Le為蒸發(fā)段長(zhǎng)度，m；Qe為蒸發(fā)段熱功率，Wt。

當(dāng)工質(zhì)溫度上升至轉(zhuǎn)變溫度Ttr時(shí)，蒸發(fā)段蒸汽區(qū)進(jìn)入局部連續(xù)流階段，蒸汽向冷凝端擴(kuò)散并在蒸發(fā)段外冷凝，使絕熱段、冷凝段熱管與工質(zhì)溫度上升。而由于熱管壁較薄，軸向熱傳導(dǎo)的效率遠(yuǎn)不如蒸汽傳熱。在熱管啟動(dòng)試驗(yàn)中發(fā)現(xiàn)，隨著熱量持續(xù)供給，蒸發(fā)段平均壁溫Taw上升至某一中間值后將保持幾乎不變，期間一條有陡峭溫度落差的蒸汽界面自蒸發(fā)段（熱區(qū)）向冷凝段（冷區(qū)）移動(dòng)，當(dāng)熱區(qū)充滿整個(gè)熱管后Taw才繼續(xù)開始上升，這一現(xiàn)象稱為熱管的“平面前鋒”現(xiàn)象，如圖5所示。

圖5 標(biāo)準(zhǔn)熱管“平面前鋒”啟動(dòng)模型示意圖Fig.5 Diagram of the standard heat pipe "flat front" startup model

Cao和Faghri據(jù)此提出了“平面前鋒”分析解模型［16］，由關(guān)于熱區(qū)長(zhǎng)度與溫降的兩個(gè)分析解組成。熱區(qū)長(zhǎng)度l是時(shí)間和熱區(qū)平均溫度Taw的函數(shù)，如下式所示：

式中：h*=2πrwh，Q（τ）=Q（τ）+h（Taw-Ta）（Le+La），ΔT=2（Taw-Ttr）；l為熱區(qū)長(zhǎng)度；La為絕熱段長(zhǎng)度；Taw為熱區(qū)平均壁溫；Ta為熱管初始溫度；Q為輸入熱量；τ為加熱時(shí)間；τi為熱區(qū)到達(dá)冷凝段時(shí)的初始時(shí)間；hfg為工質(zhì)汽化潛熱；R為氣體常數(shù)；Trf與Prf為工質(zhì)飽和狀態(tài)參考溫度與壓力；rv為蒸汽區(qū)半徑；μv為蒸汽動(dòng)力黏度。

正常工作時(shí)的熱管常使用Zuo和Faghri提出的熱阻網(wǎng)格模型［17］，該模型將熱管的每個(gè)傳熱過(guò)程視作一維熱阻導(dǎo)熱，而忽略了蒸汽的熱容和熱阻，如圖6所示。一維熱阻模型存在較多的局限性，比如未考慮可能遇到的傳熱極限問(wèn)題、忽略了絕熱段徑向?qū)崆覠o(wú)法考慮實(shí)際中絕熱段的壁面可能存在的熱損失、溫度場(chǎng)不夠精細(xì)等。為了改進(jìn)傳統(tǒng)熱阻網(wǎng)絡(luò)模型對(duì)于熱管堆系統(tǒng)瞬態(tài)分析的諸多不適用性，本文基于熱阻網(wǎng)絡(luò)模型的分析思想建立了連續(xù)流階段的二維傳熱網(wǎng)格模型。模型仍然將熱管壁、液環(huán)與吸液芯內(nèi)的導(dǎo)熱過(guò)程視作固體接觸導(dǎo)熱，但考慮了運(yùn)行中可能發(fā)生的傳熱極限問(wèn)題，將熱管吸液芯與蒸汽區(qū)間視作熱流邊界對(duì)蒸汽的換熱量單獨(dú)求解，并對(duì)于熱管結(jié)構(gòu)作出了更精細(xì)的控制體劃分，如圖6所示。模型中的各節(jié)點(diǎn)溫度的控制方程為二維導(dǎo)熱全隱差分方程，易于使用交替隱式差分（Alternating-direction Implicit，ADI）算法進(jìn)行求解［14］。

圖6 SBL-30熱管正常運(yùn)行時(shí)的傳熱模型 (a) 熱管一維熱阻網(wǎng)絡(luò)模型，(b) 改進(jìn)的二維熱阻網(wǎng)絡(luò)模型Fig.6 Heat transfer model of heat pipe in normal operation (a) Heat pipe one-dimensional thermal resistance network model,(b) Improved two-dimensional thermal resistance network model

2 KRUSTY系統(tǒng)建模與整體驗(yàn)證

2.1 建模研究與穩(wěn)態(tài)工況

基于斯特林技術(shù)的千瓦級(jí)熱管堆實(shí)驗(yàn)是目前國(guó)內(nèi)外唯一完成的帶核熱管堆全系統(tǒng)實(shí)驗(yàn)，該實(shí)驗(yàn)開展了28 h，包括了熱管堆啟動(dòng)、穩(wěn)態(tài)運(yùn)行、瞬態(tài)運(yùn)行與多種事故工況。KRUSTY熱管堆是一個(gè)小型的實(shí)驗(yàn)快中子堆，設(shè)計(jì)的穩(wěn)態(tài)熱功率為4 kW，輸出電功率1 kWe；使用三個(gè)全固態(tài)的U-8Mo環(huán)形燃料元件堆疊形成堆芯，中心為B4C控制棒；圓柱燃料側(cè)面挖槽嵌入8根鈉熱管，呈對(duì)稱式布置，熱管將熱量自堆芯傳遞給熱電轉(zhuǎn)換裝置，KRUSTY反應(yīng)堆的布置如圖7所示。因?yàn)榭紤]到成本限制與斯特林機(jī)可能存在制造差異而帶來(lái)誤差，實(shí)驗(yàn)僅安裝了兩個(gè)ASC-E2型小型斯特林機(jī)，另外使用了6個(gè)氮?dú)饫鋮s裝置來(lái)模擬斯特林機(jī)，氣冷模擬器的功率跟斯特林機(jī)保持一致。

圖7 KRUSTY堆芯布置圖 (a) 堆芯橫截面圖，(b) 堆芯實(shí)物圖Fig.7 Diagram of KRUSTY core layout (a) Cross-sectional view of reactor, (b) Physical view of reactor

KRUSTY反應(yīng)堆采用的是自通動(dòng)脈鈉熱管（Self-venting arterial heat pipe），該熱管插入堆芯后必須豎直放置，啟動(dòng)前熱管蒸發(fā)段底部聚積有鈉液池，受到加熱后產(chǎn)生鈉蒸氣，擴(kuò)散至冷凝段傳熱，鈉蒸氣液化后主要自內(nèi)壁面的動(dòng)脈通道由重力作用向下回流至液池形成循環(huán)，吸液芯網(wǎng)構(gòu)成的液體動(dòng)脈通道僅附著一部分內(nèi)壁面，如圖8所示［18］。盡管結(jié)構(gòu)與標(biāo)準(zhǔn)高溫?zé)峁芙Y(jié)構(gòu)有所區(qū)別，但瞬態(tài)運(yùn)行的機(jī)理仍然是一致的，在建模時(shí)可以按照體積等效幾何尺寸轉(zhuǎn)化為標(biāo)準(zhǔn)的壁面-液環(huán)-吸液芯三層模型。

圖8 自通動(dòng)脈鈉熱管結(jié)構(gòu)圖Fig.8 Structure diagram of self-passed arterial sodium heat pipe

為驗(yàn)證TAPIRS-D系統(tǒng)瞬態(tài)分析程序計(jì)算的準(zhǔn)確性，對(duì)該熱管堆進(jìn)行了較完整的建模，對(duì)多實(shí)驗(yàn)工況進(jìn)行了模擬計(jì)算，并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比［8］。系統(tǒng)程序中堆芯功率瞬變模型使用文獻(xiàn)中根據(jù)MCNP計(jì)算得到的點(diǎn)堆中子學(xué)參數(shù)［19］，如表1所示，并考慮了燃料、反射層和熱管溫度負(fù)反饋效應(yīng)，總的溫度負(fù)反饋系數(shù)為-0.2 $·K-1，該反應(yīng)堆沒有控制鼓系統(tǒng)，因此省略控制鼓模型。

表1 SAIRS-C熱管堆的關(guān)鍵參數(shù)Table 1 SAIRS-C heat pipe reactor key parameters

在堆芯建模時(shí)可將堆芯簡(jiǎn)化為控制棒-堆芯-熱管的三層徑向結(jié)構(gòu)，堆芯結(jié)構(gòu)與程序建模示意圖如圖9所示。將堆芯徑向劃分為兩層控制棒與三層燃料，燃料最外層（第5層）與熱管壁間設(shè)置有熱流邊界條件，熱管壁溫度為熱管蒸發(fā)段壁面節(jié)點(diǎn)的溫度，可根據(jù)傅里葉導(dǎo)熱公式計(jì)算熱流量。由于堆芯與反射層間存在真空隔熱層，堆芯向反射層的散熱較小，將堆芯對(duì)外熱損失計(jì)入第5層的內(nèi)熱源項(xiàng)。

圖9 KRUSTY堆芯的建模示意圖Fig.9 Diagram of KRUSTY core modeling

對(duì)熱管建模時(shí)，熱管幾何結(jié)構(gòu)參數(shù)由文獻(xiàn)中得到，蒸發(fā)段邊界條件為堆芯向熱管的傳遞熱量，冷凝段邊界條件為已知換熱系數(shù)與斯特林機(jī)熱端溫度。由于缺乏斯特林機(jī)與水冷模擬器的具體參數(shù)，根據(jù)文獻(xiàn)中測(cè)量的斯特林機(jī)熱端溫度數(shù)據(jù)作為熱管的邊界條件輸入程序。

本文開發(fā)的布雷頓熱電轉(zhuǎn)換模型是開式布雷頓系統(tǒng)，壓縮機(jī)入口與汽輪機(jī)出口并不相連通，未形成閉合管道回路，由于缺乏開式布雷頓系統(tǒng)相應(yīng)實(shí)驗(yàn)數(shù)據(jù)，為驗(yàn)證該模型，選擇了Sandia實(shí)驗(yàn)室的SBL-30閉式布雷頓系統(tǒng)的運(yùn)行實(shí)驗(yàn)進(jìn)行建模計(jì)算［19］，將閉式布雷頓系統(tǒng)中部分參數(shù)作為邊界條件輸入，近似處理為開式布雷頓系統(tǒng)。

整個(gè)KRUSTY反應(yīng)堆系統(tǒng)的主要參數(shù)如表2所示。

表2 KRUSTY反應(yīng)堆部分參數(shù)Table 2 KRUSTY reactor main parameters

對(duì)KRUSTY的建模完成后，首先選擇了28 h實(shí)驗(yàn)中7.5~8 h時(shí)的穩(wěn)態(tài)運(yùn)行工況，給定各熱管冷凝段的邊界平均熱流量進(jìn)行了穩(wěn)態(tài)計(jì)算。計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對(duì)比如表3所示，程序計(jì)算得到的穩(wěn)態(tài)堆芯裂變功率和燃料表面平均溫度均與實(shí)驗(yàn)數(shù)據(jù)符合良好，誤差小于1%，證明了TAPIRS-D程序?qū)τ跓峁芏颜w穩(wěn)態(tài)傳熱過(guò)程和堆芯功率反饋模擬的準(zhǔn)確性。

表3 穩(wěn)態(tài)計(jì)算結(jié)果驗(yàn)證Table 3 Steady state result comparison

2.2 熱管堆冷啟動(dòng)工況

KRUSTY實(shí)驗(yàn)的0~1 h時(shí)進(jìn)行了反應(yīng)堆冷啟動(dòng)實(shí)驗(yàn)，初始時(shí)反應(yīng)堆與環(huán)境達(dá)到熱平衡，斯特林機(jī)處于關(guān)閉狀態(tài)，在T=0.37 h時(shí)，通過(guò)引入反應(yīng)性使堆芯功率由0開始逐漸提升至最高3.7 kW并穩(wěn)定在3.1 kW，堆芯與熱管壁溫由初始環(huán)境溫度300 K上升至800 K左右，在T=0.92 h時(shí)熱管成功啟動(dòng)進(jìn)入穩(wěn)定工作狀態(tài)。

文獻(xiàn)［12］給出了一個(gè)靠近堆芯上部的熱管壁測(cè)溫點(diǎn)（測(cè)點(diǎn)1）和一個(gè)熱管絕熱段中部測(cè)溫點(diǎn)處8根熱管的平均壁溫值（測(cè)點(diǎn)2）。根據(jù)熱管堆啟動(dòng)初始狀態(tài)與功率變化數(shù)據(jù)，程序計(jì)算的熱管測(cè)點(diǎn)處壁溫結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比如圖10所示?？梢钥吹剑?.4～0.6 h期間，測(cè)點(diǎn)1處的溫度開始逐漸上升，這表明反應(yīng)堆的熱量已經(jīng)逐漸傳遞到與反應(yīng)堆相接觸的熱管蒸發(fā)段，在0.6 h左右蒸發(fā)段溫度超過(guò)400 ℃，鈉蒸氣開始在熱管內(nèi)運(yùn)動(dòng)，在0.67 h左右，測(cè)點(diǎn)2的溫度開始出現(xiàn)了急劇增加，表明此時(shí)鈉蒸氣運(yùn)動(dòng)到了熱管絕熱段的中部。從圖10可以看出，計(jì)算的測(cè)點(diǎn)溫度值與實(shí)驗(yàn)值的變化趨勢(shì)基本一致，表明TAPIRS-D程序的模型能夠合理地反映出熱管啟動(dòng)時(shí)的壁溫變化與前文提出的“平面前鋒”的推進(jìn)情況。

熱管堆冷啟動(dòng)計(jì)算過(guò)程中熱管壁溫的偏差主要可能來(lái)源于：1）程序中堆芯與熱管材料的熱物性擬合函數(shù)與實(shí)驗(yàn)中材料實(shí)際值存在一定差異；2）KRUSTY使用的熱管是垂直放置的“自通動(dòng)脈”式熱管，底部有鈉液池，其工質(zhì)回流受到重力作用，內(nèi)部結(jié)構(gòu)也不同于傳統(tǒng)高溫?zé)峁?，盡管啟動(dòng)機(jī)理同高溫?zé)峁茴愃?，但程序使用的熱管啟?dòng)模型計(jì)算會(huì)存在一定偏差。

2.3 負(fù)載變化瞬態(tài)工況

實(shí)驗(yàn)中T=8.0 h時(shí)將氣冷模擬器的流量下降一半，使每個(gè)模擬器的冷卻功率由290 W降低至210 W左右，以模擬熱管堆負(fù)載降低瞬態(tài)。當(dāng)初始時(shí)刻降低熱管冷卻功率時(shí)，熱管與堆芯由于熱量累計(jì)導(dǎo)致溫度上升，引起溫度負(fù)反饋效應(yīng)，使堆芯功率下降而遏制溫度上升，只造成了很小的溫度波動(dòng)，最終震蕩衰減至達(dá)到一個(gè)新的穩(wěn)態(tài)，堆芯功率由2.75 kW下降至2.05 kW左右。TAPIRS-D程序通過(guò)改變熱管邊界條件計(jì)算了該瞬態(tài)工況的堆芯功率與系統(tǒng)溫度響應(yīng)，計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對(duì)比如圖11所示。

圖11 負(fù)載降低工況驗(yàn)證結(jié)果 (a) 堆芯裂變功率瞬態(tài)變化，(b) 燃料組件溫度瞬態(tài)變化Fig.11 Verification results of load reduction condition (a) Core fission power transients, (b) Fuel assembly temperature transients

在T=9 h時(shí)，氣冷模擬器的冷卻功率調(diào)回正常值，反應(yīng)堆重新進(jìn)入穩(wěn)態(tài)工況，裂變功率回到2.75 kW。隨后，當(dāng)T=10.02 h時(shí)，氣冷模擬器的流量上升一倍，使每個(gè)模擬器的冷卻功率由295 W提高至510 W左右，以模擬負(fù)載增大瞬態(tài)。當(dāng)提高熱管冷卻功率時(shí)，熱管與堆芯的溫度下降使堆芯反應(yīng)性上升，導(dǎo)致堆芯功率上升，產(chǎn)生堆芯功率與溫度的震蕩并逐漸衰減至達(dá)到一個(gè)新的穩(wěn)態(tài)，堆芯功率由2.75 kW上升至2.05 kW左右。程序計(jì)算了該瞬態(tài)工況的堆芯功率與系統(tǒng)溫度響應(yīng)，與實(shí)驗(yàn)數(shù)據(jù)對(duì)比如圖12所示。

圖12 負(fù)載增大工況驗(yàn)證結(jié)果 (a) 堆芯裂變功率瞬態(tài)變化，(b) 燃料組件溫度瞬態(tài)變化Fig.12 Validation results for load increase conditions (a) Core fission power transients, (b) Fuel assembly temperature transients

由于熱管堆的負(fù)反饋效應(yīng)，在負(fù)載發(fā)生降低或者增大工況時(shí)堆芯功率都能夠收斂至穩(wěn)態(tài)，從上述分析可以看到，新的穩(wěn)態(tài)下TAPIRS-D程序預(yù)測(cè)得到的穩(wěn)態(tài)功率與溫度水平與實(shí)驗(yàn)數(shù)據(jù)符合較好，計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)震蕩的周期也比較一致，但在功率和溫度的振幅方面誤差較大，這可能是因?yàn)镵RUSTY建模信息不夠完整導(dǎo)致的?？偟膩?lái)說(shuō)，程序預(yù)測(cè)的系統(tǒng)參數(shù)響應(yīng)趨勢(shì)與實(shí)驗(yàn)數(shù)據(jù)一致，溫度的最大誤差小于1%，功率的最大誤差小于10%，計(jì)算的功率與溫度峰值與實(shí)驗(yàn)數(shù)據(jù)峰值時(shí)間對(duì)應(yīng)基本一致，最終預(yù)測(cè)應(yīng)達(dá)到的穩(wěn)態(tài)位置也與實(shí)驗(yàn)數(shù)據(jù)符合良好。

2.4 熱管失效事故工況

根據(jù)熱管在燃料芯塊的角度位置將燃料劃分為8個(gè)區(qū)域，實(shí)驗(yàn)中T=12.0 h時(shí)，堆芯0o方向的熱管所連接的模擬器氮?dú)饬髁勘磺袛?，以模擬單點(diǎn)故障導(dǎo)致的熱管失效事故，該處熱管的散熱功率由290 W降低至120 W。隨后在T=12.5 h時(shí)，增加了0o位置以外的模擬器和斯特林機(jī)的功率，平均散熱功率增加約50 W，模擬了故障后需要增加負(fù)載使總發(fā)電功率回到正常水平的操作。實(shí)驗(yàn)測(cè)量了與0o區(qū)域相鄰的45o區(qū)域和相對(duì)的180o區(qū)域的燃料外層溫度和堆芯裂變功率的響應(yīng)情況，程序計(jì)算了這一事故工況，結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對(duì)比如圖13所示。

圖13 熱管失效事故工況驗(yàn)證結(jié)果 (a) 堆芯裂變功率瞬態(tài)變化，(b) 燃料組件溫度瞬態(tài)變化Fig.13 Validate results of heat pipe failure accident (a) Core fission power transients, (b) Fuel assembly temperature transients

程序計(jì)算的功率與溫度響應(yīng)與實(shí)驗(yàn)數(shù)據(jù)符合較好，瞬態(tài)變化趨勢(shì)一致，溫度的最大誤差小于1%，功率的最大誤差小于4%。溫度計(jì)算的誤差主要來(lái)源于初始取值與實(shí)驗(yàn)值的偏差，這是因?yàn)閷?shí)驗(yàn)中堆芯功率存在周向不均勻分布，導(dǎo)致燃料芯塊周向溫度也存在不均勻分布的原因。但由于文獻(xiàn)中缺乏三維功率分布數(shù)據(jù)，程序無(wú)法模擬這種不均勻分布情況，故計(jì)算時(shí)燃料的初始穩(wěn)態(tài)溫度取兩位置測(cè)量值的平均值。當(dāng)0o位置熱管散熱功率降低時(shí)，該處燃料的多余熱量只能從其他位置的熱管排出，使相鄰的45o區(qū)域局部溫度上升；但同時(shí)堆芯熱量導(dǎo)出受阻導(dǎo)致了堆芯平均溫度升高，燃料的溫度反饋引入了負(fù)反應(yīng)性使堆芯整體功率降低。而由于0o區(qū)域向相對(duì)的180o區(qū)域傳遞熱量很小，故180o區(qū)域燃料的局部溫度反而因堆芯功率的下降而降低，使堆芯周向溫度分布變得更加不均勻。在T=12.5 h時(shí)，散熱功率上升導(dǎo)致堆芯溫度下降，引起了堆芯功率的上升，會(huì)使熱管失效區(qū)域的局部溫度繼續(xù)上升。程序計(jì)算結(jié)果很好地表現(xiàn)出了上述溫度分布不均勻性的現(xiàn)象。

2.5 反應(yīng)性引入事故工況

在實(shí)驗(yàn)的T=15.0~16.0 h，反應(yīng)堆調(diào)回至穩(wěn)態(tài)功率水平，T=16.0~16.05 h時(shí)，通過(guò)控制堆芯外圍徑向反射層套筒向下移動(dòng)0.5 mm，引入了5.5 $的負(fù)反應(yīng)性。使用程序從穩(wěn)態(tài)開始計(jì)算了該負(fù)反應(yīng)性引入事故工況，結(jié)果如圖14所示。

圖14 負(fù)反應(yīng)性引入事故下斯特林熱電轉(zhuǎn)換的參數(shù)響應(yīng) (a) 堆芯裂變功率瞬態(tài)變化，(b) 燃料組件溫度瞬態(tài)變化Fig.14 Parametric response of Stirling thermoelectric conversion under negative reactivity introduction accident(a) Core fission power transients, (b) Fuel assembly temperature transients

反應(yīng)堆在引入負(fù)反應(yīng)性后功率迅速下降導(dǎo)致溫度下降，溫度的負(fù)反饋效應(yīng)迅速補(bǔ)償了該負(fù)反應(yīng)性使功率回升，經(jīng)一定震蕩后功率回到正常穩(wěn)態(tài)水平，計(jì)算表明最終燃料平均溫度下降了27.7 K，這與實(shí)驗(yàn)測(cè)量的前后溫差29 K基本一致。

在實(shí)驗(yàn)的T=18.0~19.0 h，反應(yīng)堆的反射層套筒從正常位置下方1.0 mm開始提升了1.5 mm，總共引入了大約17.8 $的正反應(yīng)性，此階段堆芯功率最大上升至5 kW左右，是整個(gè)KRUSTY熱管堆實(shí)驗(yàn)中的峰值功率。正反應(yīng)性引入事故工況的計(jì)算結(jié)果如圖15所示，反應(yīng)堆在引入正反應(yīng)性后功率迅速上升，溫度升高引起負(fù)反應(yīng)性補(bǔ)償又使功率回降，震蕩衰減后達(dá)到新的穩(wěn)態(tài)。

圖15 正反應(yīng)性引入工況驗(yàn)證結(jié)果 (a) 堆芯裂變功率瞬態(tài)變化，(b) 燃料組件溫度瞬態(tài)變化Fig.15 Validate results of positive reactivity introduction(a) Core fission power transients, (b) Fuel assembly temperature transients

程序預(yù)測(cè)的功率和溫度響應(yīng)與實(shí)驗(yàn)數(shù)據(jù)的趨勢(shì)基本保持了一致，盡管變化振幅和功率峰值谷值仍然較實(shí)驗(yàn)數(shù)據(jù)有一定偏離，但對(duì)于燃料溫度的預(yù)測(cè)偏差較小，溫度的最大誤差小于2%，最大功率計(jì)算誤差約為16%，功率平均計(jì)算誤差小于10%。

本文在進(jìn)行實(shí)驗(yàn)驗(yàn)證時(shí)，無(wú)法準(zhǔn)確獲得燃料或包殼的真實(shí)物性，而是通過(guò)實(shí)驗(yàn)文獻(xiàn)對(duì)燃料或包殼的描述，預(yù)估出物性參數(shù)用于瞬態(tài)計(jì)算。因此，在瞬態(tài)計(jì)算過(guò)程中，燃料溫度變化的幅度會(huì)與真實(shí)情況有差異，同時(shí)，燃料溫度變化通過(guò)多普勒效應(yīng)，進(jìn)一步影響反應(yīng)堆的功率變化，從而導(dǎo)致反應(yīng)堆的堆芯功率和燃料溫度的震蕩幅度與實(shí)際測(cè)量存在一定偏差。

2.6 熱阱喪失事故工況

實(shí)驗(yàn)T=20 h時(shí)，將斯特林機(jī)空轉(zhuǎn)，模擬器氣體流量減至最低，使全部熱管僅靠環(huán)境中空氣的對(duì)流進(jìn)行被動(dòng)散熱，模擬了反應(yīng)堆的熱阱部分喪失事故，該工況下熱管冷凝端的平均冷卻功率由270 W下降至約130 W，程序模擬了這一事故工況，計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對(duì)比如圖16所示。熱阱冷卻功率的下降導(dǎo)致熱管與堆芯溫度快速上升，燃料的溫度負(fù)反饋效應(yīng)又使堆芯功率降低，遏制了堆芯溫度的上升并最終達(dá)到與降低后的冷卻功率相匹配的新穩(wěn)態(tài)。計(jì)算的燃料外層溫度在2 h前后上升了約4 ℃，其中，2 ℃是由于實(shí)驗(yàn)的溫度漂移現(xiàn)象導(dǎo)致。此外，由于低功率下的堆芯徑向溫度分布較高功率下更加均勻，故燃料外層溫度較內(nèi)層上升更快，這與實(shí)驗(yàn)測(cè)量結(jié)果基本一致，溫度的最大計(jì)算誤差小于1%，平均功率誤差小于5%。

圖16 熱阱喪失事故工況結(jié)果驗(yàn)證結(jié)果 (a) 堆芯裂變功率瞬態(tài)變化，(b) 燃料組件溫度瞬態(tài)變化Fig.16 Validation results for loss of heat sink accident (a) Core fission power transients, (b) Fuel assembly temperature transients

3 結(jié)語(yǔ)

本文使用自主研發(fā)的熱管堆瞬態(tài)分析程序TAPIRS-D對(duì)KRUSTY熱管堆實(shí)驗(yàn)系統(tǒng)進(jìn)行了完整的建模，并模擬了KRUSTY實(shí)驗(yàn)開展的熱管堆冷啟動(dòng)、運(yùn)行穩(wěn)態(tài)、負(fù)載變化、熱管失效事故、反應(yīng)性引入事故與熱阱喪失事故工況實(shí)驗(yàn)，得到了各工況下堆芯功率、溫度與反應(yīng)性的瞬態(tài)響應(yīng)，并與實(shí)際測(cè)量數(shù)據(jù)進(jìn)行了對(duì)比。結(jié)果表明：程序預(yù)測(cè)的系統(tǒng)響應(yīng)與實(shí)驗(yàn)數(shù)據(jù)符合較好，計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對(duì)比符合較好，各工況下程序計(jì)算的燃料溫度的最大相對(duì)誤差小于2%，整體功率的平均相對(duì)誤差小于10%，預(yù)測(cè)的堆芯功率與溫度瞬態(tài)響應(yīng)與實(shí)驗(yàn)數(shù)據(jù)趨勢(shì)符合一致。驗(yàn)證結(jié)果表明，自主研發(fā)的TAPIRS-D在堆芯功率瞬變模型、堆芯傳熱模型與熱管傳熱模型方面的準(zhǔn)確性與可靠性，可為熱管堆熱工與安全分析提供有效參考。

作者貢獻(xiàn)聲明吳攀負(fù)責(zé)程序研發(fā)、數(shù)據(jù)收集及系統(tǒng)建模；歐陽(yáng)澤宇負(fù)責(zé)程序研發(fā)、數(shù)據(jù)收集及系統(tǒng)建模；朱煜負(fù)責(zé)系統(tǒng)建模與驗(yàn)證；單建強(qiáng)提供技術(shù)材料支持并指導(dǎo)。