熱管冷卻反應堆堆芯瞬態(tài)熱力耦合研究

劉 博 劉利民 丁冠群 肖 瑤 顧漢洋

（上海交通大學核科學與工程學院 上海200240）

隨著世界各國對空間、深海探測需求的逐漸提高，相對于傳統(tǒng)動力方式，能量密度更高、適應環(huán)境能力更強的核動力電源成為了航天和深海探測任務的更優(yōu)選擇。熱管冷卻反應堆是目前引起關(guān)注的核反應堆電源之一，燃料裂變產(chǎn)生的熱量通過熱管傳導至熱電轉(zhuǎn)換裝置［1］，與常規(guī)反應堆相比結(jié)構(gòu)更簡單緊湊，省略了泵、閥等部件，熱管散熱具有非能動性和優(yōu)良啟動特性［2-3］，能夠?qū)崿F(xiàn)更高的可靠性和安全性［4］。

在高溫運行時熱管冷卻反應堆的熱膨脹效應和熱應力效應顯著，因此在堆芯設計階段必須考慮熱工、應力安全問題。Kapernick 等［5］針對SAFE（Safe Affordable Fission Engine）熱管冷卻反應堆（SAFE-100）通過有限元軟件COSMOS 和SINDA/FLUENT耦合計算獲得堆芯溫度場和熱應力。Wang 等［6］提出了小型熱管冷卻反應堆（small Heat Pipe cooled Reactor，sHPR），開發(fā)了熱管冷卻反應堆熱工水力分析程序用于分析sHPR 的穩(wěn)態(tài)和瞬態(tài)性能，得到各關(guān)鍵部件熱工參數(shù)，并對熱管冷卻反應堆的典型事故進行了分析和評價，包括功率變化、單根熱管失效、冷卻通道堵塞等情況。張文文等［7］針對新型空間熱管反應堆采用計算流體力學軟件FLUENT 對堆芯進行了穩(wěn)態(tài)熱工安全分析，對控制轉(zhuǎn)鼓7 種不同角度下的正常工況以及單根熱管失效的事故工況進行計算分析，得到最熱通道各層材料的溫度分布。馬譽高等［8］提出了一種考慮熱管堆固體堆芯顯著膨脹的幾何更新和反應性反饋方法，并構(gòu)建了基于動態(tài)幾何的中子物理/熱工/力學三場核熱力耦合分析程序，基于此方法對MegaPower 熱管堆進行了核熱力耦合分析。

雖然關(guān)于熱管冷卻反應堆熱工、力學特性的問題已有很多研究，但大多數(shù)研究沒有考慮核熱力的耦合效應，也很少有堆芯事故瞬態(tài)變化情況的研究?；谝陨蠁栴}，本文基于有限元分析軟件FLUENT以及MECHANICAL針對堆芯的典型發(fā)熱單元分別進行核熱力耦合分析，采用點堆中子動力學模型進行堆芯功率、溫度場的瞬態(tài)計算；并最終基于溫度場和力學場分布特征開展堆芯設計以及事故工況下的安全評估。

1 耦合數(shù)值方法

1.1 核熱耦合模型

1.1.1 點堆中子動力學模型

對NUSTER（Nuclear Silence ThermoElectric Reactor）反應堆進行瞬態(tài)核熱耦合分析時，反應堆燃料裝載量較少，堆芯狀態(tài)位于臨界附近且中子密度隨時間變化較緩慢，因此，可采用考慮6組緩發(fā)中子的點堆動態(tài)方程來求解瞬時燃料棒裂變功率［9］。點堆動力學模型描述了反應堆中子密度與反應性之間的瞬時變化關(guān)系，該模型忽略空間效應，計算簡便、響應快速。點堆動力學模型見式（1）、（2）：

式中：P(t)為裂變功率，W；ρ(t)為隨時間變化的反應性，$；β為緩發(fā)中子份額；Λ為平均中子代時間，s；λ為衰變常數(shù)；C(t)為緩發(fā)中子先驅(qū)核濃度；下標i代表第i組緩發(fā)中子。

使用全隱式一階泰勒多項式積分方法對點堆動力學方程進行求解，既可以解決剛性問題，又具有精度高、適用性強的特點［10］。因此，本文采用全隱式一階泰勒多項式積分法，將顯式數(shù)值求解后的點堆中子動力學模型基本方程編入FLUENT 軟件的用戶自定義函數(shù)（User Defined Function，UDF）中，以實現(xiàn)在FLUENT軟件內(nèi)的核熱耦合迭代計算。

1.1.2 反應性反饋模型

本文研究對象物理模型考慮的反應性反饋主要包括燃料反應性反饋以及導熱基體反應性反饋，計算過程中的系統(tǒng)反應性變化如式（3）、（4）所示：

式中：ρ(t)為總反應性，$；ρ0為初始反應性；ρex(t)為外部引入反應性，$；ρi(t)為堆內(nèi)各材料反饋反應性，$，包括燃料多普勒效應反應性反饋ρU和基體膨脹效應反應性反饋ρMo。

UO2燃料芯塊中，燃料溫度效應主要是由燃料核共振吸收的多普勒效應引起的。溫度升高使燃料中有效共振吸收增加，有效增殖系數(shù)下降，產(chǎn)生負的反應性反饋，如式（5）所示：

式中：αDoppler為多普勒溫度系數(shù)，$·K-1；TF，avg(t)為t時刻燃料棒平均溫度，K；TF，avg(0)為正常穩(wěn)態(tài)時的燃料棒平均溫度，K。

導熱基體由于溫度變化會產(chǎn)生材料密度變化、熱膨脹造成的尺寸變化等效應，這些效應產(chǎn)生的基體反應性反饋由式（6）給出：

式中：αMo為基體溫度反應性系數(shù)，$·K-1；TMo，avg(t)為t時刻基體的平均溫度，K；TMo，avg(0)為正常穩(wěn)態(tài)時的基體平均溫度，K。

1.2 熱力耦合模型

熱管堆區(qū)別于其他傳統(tǒng)堆型最顯著的特征是其堆芯的固態(tài)屬性，因此熱膨脹效應和熱應力問題是其研究的重點［1］。熱應力問題的研究主要在于確定溫度場以及由溫度場確定應力應變。物體內(nèi)部的溫度差ΔT(x，y，z)將引起αTΔT(x，y，z)的熱膨脹，則該彈性物體的物理方程變成了式（7）所示：

式中：σi為應力分量，Pa；εi為形變分量，m·m-1；αT為熱膨脹系數(shù)，K-1；E為彈性模量，Pa；G為切變模量，Pa；μ為泊松系數(shù)。

2 計算模型及參數(shù)

2.1 幾何模型及邊界條件

本文研究對象為一種典型熱管冷卻反應堆——靜默式海洋熱管冷卻反應堆［11］，堆芯主要包括導熱基體、燃料棒、熱管、軸向反射層、滑動反射層、控制棒組件等［12］。NUSTER堆芯可劃分為109個由燃料棒、導熱基體以及中心熱管構(gòu)成的發(fā)熱單元［13］，發(fā)熱單元于堆芯中的選取方式及發(fā)熱單元示意圖如圖1所示。由于開展全堆芯模擬計算量巨大，故建立發(fā)熱單元的三維模型如圖2 所示作為研究對象；為了完整模擬熱管區(qū)域的傳熱，對熱管蒸發(fā)段以及堆芯外部分的熱管絕熱段、冷凝段進行建模，模擬堆芯到冷源的完整傳熱過程。NUSTER堆芯燃料區(qū)內(nèi)部分結(jié)構(gòu)尺寸參數(shù)以及組件材料選擇如表1所示。

表1 燃料區(qū)內(nèi)結(jié)構(gòu)參數(shù)Table 1 Design parameters of fuel area

圖1 發(fā)熱單元示意圖Fig.1 Schematic diagram of heat unit

圖2 發(fā)熱單元三維模型Fig.2 Schematic diagram of three-dimensional model of heat unit

根據(jù)實際情況和物理意義，邊界條件設置如下：熱分析中邊界條件為熱管冷凝段壁面為定壁溫條件。力學分析中邊界條件設置為側(cè)面為對稱邊界條件，兩個頂面和底面為軸向位移約束邊界條件，如圖3所示。

圖3 力學邊界條件Fig.3 Mechanical boundary condition

2.2 物性參數(shù)

在熱分析中涉及的各部件材料熱物性相關(guān)參數(shù)［14-16］如表2 所示，物性在FLUENT 軟件中由用戶自定義程序UDF編寫給出。

表2 材料熱物性參數(shù)Table 2 Thermophysical properties of materials

其中熱管蒸汽區(qū)導熱系數(shù)［15-16］如式（8）所示：

式中：k為蒸汽導熱系數(shù)，W·（m·K）-1；R為蒸汽區(qū)熱阻，K·W-1；Leq為熱管等效長度，m；μv為蒸汽動力黏度，Pa·s；Tv為蒸汽溫度，K；dv為蒸汽區(qū)直徑，m；ρv為蒸汽密度，kg·m-3；hfg為汽化潛熱，J·kg-1。

在MECHANICAL軟件中設置各固體材料的力學物性參數(shù)，UO2熱膨脹關(guān)系與塑性區(qū)的應變應力方程［17］分別如式（9）、（10）所示，其余力學物性參數(shù)如表3所示。

表3 材料力學物性參數(shù)Table 3 Mechanical properties of materials

式中：T為溫度，K，范圍為300~2 873 K；Δl為相對于 300 K時的線性熱膨脹。

式中：σ為塑性變形階段的應力，MPa；εp為塑性變形 階段的應變，m·m-1。

瞬態(tài)計算中，考慮核反應堆中子動力學和多普勒反應性反饋、導熱基體溫度反應性反饋等，以及通過堆芯物理計算得到的NUSTER 的堆芯反應性反饋系數(shù)表達式如表4所示。

表4 NUSTER中子動力學參數(shù)Table 4 Neutron dynamic parameters of NUSTER

2.3 網(wǎng)格劃分

熱分析采用六面體網(wǎng)格劃分，軸向網(wǎng)格拉伸，力分析采用六面體與四面體結(jié)合的網(wǎng)格劃分；發(fā)熱單元網(wǎng)格劃分示意圖分別如圖4（a）、（b）所示。兩種網(wǎng)格劃分均滿足網(wǎng)格無關(guān)性，網(wǎng)格敏感性計算結(jié)果分別如圖5（a）、（b）所示，據(jù)此熱分析網(wǎng)格數(shù)量約為200萬，力學分析網(wǎng)格數(shù)量約為28萬。

圖4 網(wǎng)格劃分示意圖 (a)熱分析網(wǎng)格，(b)應力分析網(wǎng)格Fig.4 Schematic of mesh (a)Mesh of heat transfer analysis,(b)Mesh of mechanical analysis

圖5 網(wǎng)格敏感性結(jié)果Fig.5 Results of mesh sensitivity

3 熱力耦合結(jié)果

NUSTER 系統(tǒng)采用4 組滑動反射層和4 個安全控制棒作為初級反應性控制系統(tǒng)，當控制系統(tǒng)發(fā)生故障或者外界因素干擾時，控制棒可能會意外抽出（提棒事故）或失控彈出（彈棒事故）而引入正反應性。反應性引入事故是造成瞬態(tài)超功率事故的典型事故之一，因此有必要研究熱管冷卻反應堆NUSTER 在正反應性引入事故情況下的瞬態(tài)特性，對堆芯安全性進行驗證。NUSTER堆芯中控制棒的價值共約800×10-5，假定彈棒事故發(fā)生時在5 s 內(nèi)線性引入400×10-5反應性，無外界擾動。針對典型穩(wěn)態(tài)工況以及上述反應性引入事故工況，進行熱力耦合分析。

3.1 穩(wěn)態(tài)工況結(jié)果

反應性引入事故發(fā)生前穩(wěn)態(tài)工況的發(fā)熱單元溫度場結(jié)果如圖6所示，溫度峰值為1 516 K。此時發(fā)熱單元和各組件的等效應力場如圖7 所示，圖7（a）為發(fā)熱單元的總體等效應力云圖剖面圖，峰值為196 MPa 位于包殼；圖7（b）為包殼組件的等效應力云圖，顯示等效應力峰值位于包殼內(nèi)邊緣尖角。圖7（c）為導熱基體的等效應力分布情況，峰值同樣位于基體外表面的棱角處，發(fā)生了應力集中，峰值為162 MPa，其余位置等效應力為120 MPa 左右。圖7（d）展示了熱管壁的等效應力分布情況，峰值為177 MPa，位于熱管內(nèi)壁面。

圖6 穩(wěn)態(tài)工況溫度云圖Fig.6 Temperature contour under steady condition

圖7 穩(wěn)態(tài)工況等效應力場云圖 (a)發(fā)熱單元，(b)包殼，(c)導熱基體，(d)熱管壁Fig.7 Equivalent stress contour under steady condition (a)Heat unit,(b)Claddings,(c)Matrix,(d)Heat pipe wall

3.2 事故工況結(jié)果

堆芯反應性和功率以及主要組件溫度峰值的變化趨勢分別如圖8所示。200 s時開始引入外部反應性，燃料棒功率隨著總反應性的增加而增大，芯塊溫度也逐漸升高。各組件溫度升高后，在燃料負反饋效應和基體負反饋效應引入的負反應性作用下，總反應性約在220 s 達到峰值后逐漸降低。由于堆芯導熱相對于反應性和功率的瞬變更加緩慢，有一定的滯后性，反應性會繼續(xù)下降至負值；因此各組件溫度會緩慢升高至約250 s達到峰值，之后溫度開始小幅度降低。最終約在350 s時，反應性、功率、組件溫度等參數(shù)達到新的穩(wěn)定狀態(tài)。引入正反應性事故下八分之一堆芯最終的穩(wěn)定功率為194 kW，約為初始穩(wěn)定狀態(tài)的1.5倍。

圖8 反應性(a)和功率(b)瞬時變化Fig.8 Transient variation of core reactivity(a)and core power(b)

圖9給出了從穩(wěn)態(tài)工況到事故發(fā)生后主要組件的溫度峰值變化情況。結(jié)果表明：事故過程中燃料芯塊溫度峰值1 700 K，比初始穩(wěn)態(tài)升高了220 K；包殼溫度峰值1 415 K，比初始穩(wěn)態(tài)升高了109 K；熱管壁1 364 K，比初始穩(wěn)態(tài)升高了90 K。將各組件溫度峰值與各材料的熔點進行比較，結(jié)果匯總于表5，結(jié)果表明各組件距離熔點均有200 K 以上的安全裕度，事故工況下的堆芯熱工安全得到驗證。

表5 事故工況下各組件溫度峰值Table 5 Peak temperature value of each assembly under accident condition

由圖9 可得，260 s 時發(fā)熱單元各組件出現(xiàn)溫度峰值，基于熱應力的保守性計算選取此時刻展示溫度與熱應力的分布情況，溫度云圖以及等效應力云圖分別如圖10、11所示。通過比較不同組件在穩(wěn)態(tài)工況和引入事故工況的等效應力峰值，得到反應性引入事故對應力分布的影響并進行安全評估。此時單元溫度峰值為1 700 K，等效應力峰值為334 MPa，位于導熱基體。圖11（a）為發(fā)熱單元等效應力剖面圖，云圖表明基體與熱管壁的等效應力比其他組件高；圖11（b）為包殼組件等效應力云圖，峰值為327 MPa；圖11（c）為導熱基體等效應力分布情況，應力峰值為334 MPa；圖11（d）為熱管壁等效應力云圖，峰值為300 MPa。對比圖7、11 的應力結(jié)果發(fā)現(xiàn)，事故發(fā)生后基體等效應力峰值最多增大172 MPa，包殼等效應力最多增大131 MPa，熱管壁峰值增大了123 MPa。

圖9 組件溫度峰值變化Fig.9 Variation of assembly peak temperature

圖10 引入反應性事故工況（t=260 s時）溫度云圖Fig.10 Temperature contour under reactivity insertion accident condition(t=260 s)

圖11 反應性引入事故（t=260 s時）等效應力場云圖 (a)發(fā)熱單元，(b)包殼，(c)導熱基體，(d)熱管壁Fig.11 Equivalent stress contour under reactivity insertion accident condition(t=260 s)(a)Heat unit,(b)Claddings,(c)Matrix,(d)Heat pipe wall

為了更細致地研究事故對不同組件等效應力分布的影響，分別沿包殼和熱管壁選擇90°、360°的周向路徑如圖12（a）、（b）所示，另外選擇徑向截面上燃料芯塊到熱管壁的徑向路徑如圖12（c）所示。

圖12 周向路徑選取示意圖 (a)包殼周向路徑，(b)熱管壁周向路徑，(c)單元徑向路徑Fig.12 Schematic diagram of Circumferential paths selection (a)Circumferential path along cladding,(b)Circumferential path along heat pipe wall,(c)Radial path in unit

分別作這三條路徑上兩種工況下的等效應力分布曲線，對比結(jié)果分別如圖13（a）、（b）、（c）所示。圖13（a）結(jié)果表明，引入反應性事故會造成包殼上等效應力的增大，在30°、60°位置等效應力差值最大，約為200 MPa。熱管壁上周向路徑的等效應力對比結(jié)果如圖13（b）所示，事故工況下的熱管壁等效應力明顯大于穩(wěn)態(tài)工況下的，在kπ/2（k=0，1，2，3）處等效應力差值最大，約203 MPa。從燃料芯塊到熱管壁徑向路徑上的等效應力分布如圖13（c）所示，結(jié)果表明反應性引入事故下沿徑向路徑上各組件的等效應力均大于穩(wěn)態(tài)工況下，燃料芯塊內(nèi)區(qū)等效應力差值約205 MPa，熱管壁差值最大約156 MPa。

圖13 反應性引入事故與穩(wěn)態(tài)工況等效應力比較(a)沿包殼周向路徑曲線，(b)沿熱管壁周向路徑曲線，(c)沿發(fā)熱單元徑向路徑曲線Fig.13 Comparison of equivalent stress between reactivity insertion accident and steady condition (a)Curve of circumferential path along cladding,(b)Curve of circumferential path along heat pipe wall,(c)Curve of radial path in heat unit

基于上述溫度、應力結(jié)果，進行事故工況下的堆芯安全評估。各組件溫度峰值變化情況已由圖10給出，表5 也表明各組件溫度均遠小于相應材料的熔點，事故下堆芯熱工安全得到驗證。依據(jù)彈塑性力學強度理論中的von-Mises屈服條件［21］，組件應力峰值與材料屈服強度對比結(jié)果總結(jié)于表6 中，關(guān)鍵組件等效應力峰值均遠小于相應材料的屈服強度，符合von-Mises屈服條件；并且組件應力距限值均有150 MPa 以上的安全裕度，因此組件應力方面滿足強度要求。

4 結(jié)語

基于有限元分析軟件ANSYS平臺，結(jié)合點堆中子動力學模型，針對兩種工況——穩(wěn)態(tài)工況與反應性引入事故工況，對靜默式海洋熱管冷卻反應堆NUSTER堆芯中的典型發(fā)熱單元開展了核熱耦合以及熱力耦合計算，給出了發(fā)熱單元溫度場、應力場以及反應性、功率的瞬態(tài)變化結(jié)果等。結(jié)果表明：事故工況發(fā)生后各組件的溫度、等效應力均有明顯增大，對事故工況下的發(fā)熱單元進行溫度、應力安全評估，對比結(jié)果表明各組件均滿足熱工、應力的限制要求；并且事故下各組件溫度有200 K 以上的安全裕度，等效應力距離安全限值有150 MPa 以上的安全裕度，證明了NUSTER堆芯設計在穩(wěn)態(tài)和事故工況下的可靠性和安全性。

作者貢獻聲明劉博：實施研究計算與數(shù)據(jù)處理，起草文章；劉利民：設計研究方案與指導；丁冠群：提供軟件技術(shù)指導；肖瑤：提供總體技術(shù)指導；顧漢洋：獲取研究經(jīng)費，提供平臺支持。