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

（上海交通大學(xué)核科學(xué)與工程學(xué)院 上海200240）

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

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

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

1 耦合數(shù)值方法

1.1 核熱耦合模型

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

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

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

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

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

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

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

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

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

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

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

1.2 熱力耦合模型

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

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

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

2.1 幾何模型及邊界條件

本文研究對象為一種典型熱管冷卻反應(yīng)堆——靜默式海洋熱管冷卻反應(yīng)堆［11］，堆芯主要包括導(dǎo)熱基體、燃料棒、熱管、軸向反射層、滑動反射層、控制棒組件等［12］。NUSTER堆芯可劃分為109個由燃料棒、導(dǎo)熱基體以及中心熱管構(gòu)成的發(fā)熱單元［13］，發(fā)熱單元于堆芯中的選取方式及發(fā)熱單元示意圖如圖1所示。由于開展全堆芯模擬計算量巨大，故建立發(fā)熱單元的三維模型如圖2 所示作為研究對象；為了完整模擬熱管區(qū)域的傳熱，對熱管蒸發(fā)段以及堆芯外部分的熱管絕熱段、冷凝段進(jìn)行建模，模擬堆芯到冷源的完整傳熱過程。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ù)實(shí)際情況和物理意義，邊界條件設(shè)置如下：熱分析中邊界條件為熱管冷凝段壁面為定壁溫條件。力學(xué)分析中邊界條件設(shè)置為側(cè)面為對稱邊界條件，兩個頂面和底面為軸向位移約束邊界條件，如圖3所示。

圖3 力學(xué)邊界條件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ū)導(dǎo)熱系數(shù)［15-16］如式（8）所示：

式中：k為蒸汽導(dǎo)熱系數(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è)置各固體材料的力學(xué)物性參數(shù)，UO2熱膨脹關(guān)系與塑性區(qū)的應(yīng)變應(yīng)力方程［17］分別如式（9）、（10）所示，其余力學(xué)物性參數(shù)如表3所示。

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

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

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

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

表4 NUSTER中子動力學(xué)參數(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萬，力學(xué)分析網(wǎng)格數(shù)量約為28萬。

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

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

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

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

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

3.2 事故工況結(jié)果

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

圖8 反應(yīng)性(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。將各組件溫度峰值與各材料的熔點(diǎn)進(jìn)行比較，結(jié)果匯總于表5，結(jié)果表明各組件距離熔點(diǎn)均有200 K 以上的安全裕度，事故工況下的堆芯熱工安全得到驗(yàn)證。

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

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

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

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

圖11 反應(yīng)性引入事故（t=260 s時）等效應(yīng)力場云圖 (a)發(fā)熱單元，(b)包殼，(c)導(dǎo)熱基體，(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

為了更細(xì)致地研究事故對不同組件等效應(yīng)力分布的影響，分別沿包殼和熱管壁選擇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

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

圖13 反應(yīng)性引入事故與穩(wěn)態(tài)工況等效應(yīng)力比較(a)沿包殼周向路徑曲線，(b)沿?zé)峁鼙谥芟蚵窂角€，(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

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

4 結(jié)語

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

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