熱管冷卻反應堆固態堆芯熱力耦合性能分析及結構優化

岳明楷，金 浩，劉桐蔚，李文杰，黃永忠，惲 迪，周進雄,*

(1.西安交通大學 航天航空學院 機械振動與強度國家重點實驗室，陜西 西安 710049；2.中國核動力研究設計院，四川 成都 610065；3.西安交通大學 能源與動力學院，陜西 西安 710049)

核能以其持久、經濟、安全和清潔等優勢被人們認為是當前最具開發價值和發展潛力的新型能源，其持續發展不僅能保證國民經濟的平穩發展，還對國家安全有重要的戰略意義[1]。

熱管冷卻反應堆一回路系統不采用傳統壓水堆、液態金屬冷卻反應堆等使用冷卻劑回路的布置方式，而采用熱管組件將堆芯產生的熱量傳遞到堆外的能量轉換系統。20世紀60 年代，LANL( Los Alamos National Laboratory)提出了采用高效的熱管作為導熱元件的新型空間核反應堆設計概念[2]。LANL的研究者們進一步研究設計，逐漸形成燃料棒、熱管相互間隔排布的熱管冷卻反應堆結構設計理念[3-4]。目前，已開展的熱管冷卻反應堆方案與關鍵技術研究已經證實了熱管冷卻反應堆具有固有安全特性、成熟技術基礎的可借鑒性以及較小的比質量等特性[5-6]。2018年3月，Kilopower反應堆概念的地面原型堆項目KRUSTY成功完成了有核實驗，標志著世界上第1個熱管堆的成功運行[7]。熱管堆便于小型化，核能具有可離線性、高功率密度、高可靠性、運行時間長，未來在偏遠地區供能供熱、科考供能、空間核動力等領域具有應用潛力[8-9]。

通過對熱管堆的燃料性能進行分析，可促進對熱管冷卻反應堆特性的了解，進一步可幫助我們對固態堆芯設計參數進行優化和改進，從而提高熱管堆的優勢同時盡可能地減弱其熱電轉換效率低的缺點。

在熱管堆工作壽期內固態堆芯長時間處于高溫、輻照環境工況中，涉及間隙傳熱、材料蠕變等復雜的熱力耦合行為。本文介紹Ross和Stoute經典間隙傳熱模型并對固態堆芯所用材料在高溫、輻照環境中的蠕變、腫脹等力學行為進行簡要說明。綜合考慮間隙傳熱模型及極端工況下材料的力學行為，基于ABAQUS結合GAPCON、CREEP子程序的編制進行二次開發，建立固態堆芯等效二維截面熱力耦合計算模型，并對5 a壽期內的工作性能進行數值模擬?；趨祷Ｓ邢拊椒▽虘B堆芯尺寸進行參數優化，通過深度神經網絡建立堆芯最高溫度、最大應力兩個優化目標與設計參數的代理模型。基于代理模型，使用NSGA-Ⅱ算法獲得pareto前沿解集，以獲取兩個優化目標權重相等情況下的最優堆芯設計。

1 固態堆芯燃料組件性能分析

熱管反應堆固態堆芯燃料組件由堆芯基體和按一定規律間隔排布在其中的燃料棒、熱管構成。其中1種典型的固態堆芯排布方式如圖1所示，固態堆芯由外圍中子反射層、反應性控制鼓、堆芯基體、熱管、燃料棒等部分組成。燃料棒由柱狀裂變燃料芯塊和合金包殼封裝組成。熱管是1963 年LANL開發的傳熱元件，通過管內工質的相變過程的高潛熱特性和毛細作用實現熱量從蒸發段到冷凝段的傳遞，是一種高效的被動傳熱方式。對于如圖1所示的熱管冷卻反應堆，燃料棒、熱管以蜂窩狀六邊形排布方式間隔排布在堆芯基體中，這樣的排布方式能最大程度上使燃料棒均勻地臨近熱管，使燃料棒釋放熱量高效地傳遞給熱管。熱管蒸發段在堆芯內與核燃料棒一起排布在固態堆芯基體中，冷凝段布置在熱電轉換系統中，反應堆工作時，核燃料內的裂變能量傳遞給金屬熱管的蒸發段，熱管將熱量從堆芯傳導至能量轉換系統轉換為電能或將殘余的熱量排出至最終熱阱[10-11]。

圖1 典型熱管冷卻反應堆結構Fig.1 Typical form of heat-pipe-cooled reactor structure

熱管冷卻反應堆固態堆芯工作狀態涉及傳熱、力學、核裂變反應物理等諸多領域的問題，各因素耦合相互影響，對這些問題的仿真和分析將有助于更好地設計固態堆芯。

1.1 含氣隙結構的傳熱

傳熱問題在核工程中普遍存在，是固態堆芯重要研究問題。在熱管冷卻反應堆固態堆芯中，核燃料棒中的裂變燃料芯塊發生核反應，產生的熱量依次經過燃料包殼、堆芯基體傳遞到熱管，經由熱管內傳熱介質排出堆芯。整個傳熱路徑經過各零件，零件裝配所致的間隙對傳熱過程有明顯的影響。傳熱路徑示意圖如圖2所示。

圖2 固態堆芯內傳熱示意圖Fig.2 Schematic diagram of heat transfer in solid core

發生在固體內部的傳熱問題由傅里葉熱傳導方程控制：

(1)

式中：k為熱導率；T為溫度；Q為熱源功率密度；ρ為密度；c為比熱容；t為時間。

間隙是傳熱路徑的重要組成部分，Ross-Stoute提出的間隙傳熱機制認為，透過間隙表面傳遞的熱量由各表面的溫度和等效間隙傳熱系數決定，間隙傳熱模型[12]如下：

qg=hgap(T1-T2)

(2)

式中：qg為間隙熱通量；hgap為間隙傳熱系數；T1、T2為各間隙表面溫度。其中，hgap由間隙內的傳熱機制決定，主要由三部分構成，包括間隙內氣體的導熱、間隙表面之間的輻射換熱以及可能發生的間隙閉合部分的接觸導熱，可由下式表示：

hgap=hg+hs+hr

(3)

式中：hg為間隙氣體傳熱系數；hs為間隙表面接觸傳熱系數；hr為間隙表面輻射換熱系數。

1) 間隙氣體傳熱模型

hg與氣隙的大小以及兩表面的粗糙度有關，間隙表面附近氣體溫度并不等于固體表面溫度，這種溫度不連續的現象被稱為溫度跳躍。間隙表面的粗糙度也是影響間隙氣體導熱的重要因素。

氣體傳熱系數計算模型[12]如下：

(4)

式中：kg為氣體熱導率；d為氣隙寬度；Cr為粗糙度系數；R1、R2分別為間隙表面粗糙度；g1、g2為間隙表面跳躍距離。

跳躍距離的計算由Kennard模型[13]給出：

(5)

式中：amix為混合氣體組分調節系數；Tg為混合氣體平均溫度；p為混合氣體壓力；n為混合氣體中所包含的氣體分子種類；Xi為第i種氣體的分子份額；Mi為第i種氣體的相對分子質量。

在間隙氣體導熱模型中，將氣體假設為不會流動的固體進行分析，忽略在加熱情況下可能發生流動而產生對流換熱的作用。氣體的熱導率隨著溫度、氣體成分和反應堆操作條件的變化而變化。例如對于由He、Xe、Kr等組成的混合氣體，使用下式計算熱導率：

kg=(kHe)x(kXe)y(kKr)z

(6)

式中，x、y、z為各氣體組分的摩爾分數。

2) 間隙接觸導熱模型

當間隙表面局部可能發生接觸時，接觸區域對間隙傳熱起到了重要作用。由于接觸表面并非理想光滑的，有一定的表面粗糙度，造成實際上只有接觸區域的部分點實際擠壓接觸在一起，而且還存在一定的氣隙。因此現有的接觸導熱模型均認為導熱系數與接觸壓有關[12]：

(7)

3) 輻射換熱模型

在氣隙換熱過程中，溫度較高的情況下，輻射換熱也是間隙換熱的重要組成部分。依據波爾茲曼輻射換熱定律有[12]：

(8)

(9)

式中：σ為Stephan-Boltzman常數；Fe為與兩表面發射率相關的輻射函數；e1、e2為間隙表面輻射率。

1.2 堆芯材料的蠕變

對于反應堆部件材料，由于長期處于高溫、高輻照的工況中，其形變成因包含了彈性變形、熱膨脹、蠕變和可能存在裂變產物引起腫脹等多種因素[14-19]，材料總應變為：

εtotlp=εep+εtp+εcp+εsw

(10)

式中：εtotlp為總應變；εep為彈性應變；εtp為熱膨脹應變；εcp為蠕變應變；εsw為腫脹應變。

對于UO2燃料芯塊，εcp和εsw由Allison給出了經驗模型[14]。芯塊蠕變模型為：

(11)

燃料芯塊的腫脹行為主要由固體裂變產物所導致的體積生長和裂變氣體所在燃料芯塊內部產生的氣泡導致的體積膨大構成，芯塊腫脹模型為：

(12)

(13)

e-0.016 2(2 800-T)-8.0×10-27F

(14)

對于并不發生裂變反應的包殼、基體等不銹鋼材料，不考慮腫脹，僅考慮蠕變[15]，包殼蠕變模型為：

(15)

1.3 熱力耦合問題

由于核反應堆工作環境處于高溫、高輻照工況下，材料屬性與溫度、輻照等因素相關。堆芯因溫度變化引起結構形變也使間隙特征改變，間隙熱導率對間隙屬性敏感，從而影響堆芯傳熱行為。因此對堆芯性能進行分析時，必須考慮多物理場相互耦合下溫度、變形、輻照等多種因素的相互作用，在ABAQUS中熱力耦合分析步進行逐步迭代計算。

在固態堆芯設計過程中，非線性是建模過程的難點。材料的腫脹、蠕變等機理導致材料非線性。間隙尺寸相對較小，隨著壽期的變化，部分間隙處可能閉合，導致接觸非線性。

2 固態堆芯熱力耦合分析有限元模型設置

常規的奧氏體不銹鋼和鎳基合金都是較為合適的包殼結構材料，同時考慮反應堆整體的經濟性，采用商業上技術比較成熟的316不銹鋼作為包殼材料，同時在堆芯材料選擇采用二氧化鈾作為堆芯材料。固態堆芯由燃料棒、熱管和基體組成，基體整體形狀為六棱柱，其中有排布燃料棒和熱管的通道，燃料棒和熱管均為細長的圓柱體，固態堆芯結構示意圖如圖3所示，其截面尺寸遠小于軸向方向尺寸，其幾何特點為長細比(軸向方向尺寸與截面尺寸的比值)大、具有對稱性。間隙尺寸至關重要，導致全尺寸的三維模型計算將耗費大量的時間成本，不適用于參數化生成大量數據樣本。

圖3 固態堆芯三維模型及橫截面示意圖Fig.3 3D model and cross-section diagram of solid core

出于計算成本的考慮，結合上述特點，對堆芯的橫截面進行廣義平面應變假設，并利用堆芯截面的對稱性，選取代表性柵元代表堆芯，如圖4所示。基于ABAQUS在堆芯設計尺寸允許范圍內選取1組參數采用熱力耦合分析步進行仿真計算。尺寸參數列于表1。計算堆芯在正常工況下5 a壽期內性能變化。在間隙等細節處布置了尺寸為0.4 mm的局部細化的網格。堆芯外表面設為絕熱邊界條件，采用簡化的恒溫熱管模型，燃料芯塊的熱功率為107W/m3。通過編寫子程序GAPCON加入了上述間隙傳熱模型，實現間隙傳熱行為的數值仿真，通過子程序CREEP加入了上述蠕變、腫脹模型，實現對材料蠕變和腫脹特性的模擬。計算參數列于表2[20-26]。

a——代表性柵元選取示意圖；b——設計尺寸變量示意圖圖4 計算模型示意圖Fig.4 Schematic diagram of model

表1 固態堆芯尺寸參數Table 1 Design size parameter of solid core

3 有限元結果分析

對于上述尺寸的固態堆芯結構，進行熱力耦合分析。隨著堆芯工作時間的增長，燃料芯塊的蠕變腫脹逐漸加劇，各部件體積逐漸增大，因此，堆芯內部狀態隨著壽期的增長而逐漸加劇。通過有限元數值模擬獲取其在正常工況下5 a壽期時的溫度場分布和面內應力分布，如圖5所示。

a——溫度場仿真結果；b——面內應力仿真結果圖5 壽期5 a時堆芯截面熱力耦合仿真結果Fig.5 Simulation results of temperature and stress at operation period of 5 years

由圖5a可知，堆芯整體溫度分布規律顯著，燃料芯塊中心溫度較高，呈現沿著燃料芯塊徑向環狀下降的特點，除了堆芯最外圍的燃料芯塊，其余部分的分布規律相似，幅值相近，這是由于熱管和燃料棒的分布具有周期性。堆芯最外圍相鄰的燃料芯塊的溫度顯著高于其他位置的燃料芯塊，這是由于該位置處附近的熱管平均數量較少。

由圖5b可看到堆芯面內應力的分布情況，據此能對芯塊開裂進行部分解釋。燃料芯塊中心受到向內擠壓，出現明顯的壓應力區域，而燃料芯塊外圍受拉，出現環向拉應力區域，這種應力特點是造成柱狀的燃料芯塊徑向開裂的重要原因。形成這樣應力分布的主要原因是，燃料芯塊產生的熱量要從外表面通過氣隙向外傳遞，從而導致燃料芯塊內部溫度高，膨脹較多，內外相互擠壓牽拉。除了燃料芯塊的熱膨脹，燃料芯塊的腫脹、蠕變特性導致燃料體積隨著時間逐漸增長，產生的蠕變應變、腫脹應變，導致燃料半徑逐漸增大，燃料芯塊與包殼之間的間隙逐步減小，逐漸縮小間隙，間隙傳熱系數增加，有利于間隙傳熱。不過此設計參數下的堆芯在設計的5 a壽期內，芯塊和包殼的間隙并未閉合，因而未產生接觸壓力，燃料包殼、熱管包殼和堆芯基體的截面內應力相較于燃料芯塊并不顯著。

表2 固態堆芯材料及計算參數Table 2 Material and calculation parameters of solid core

4 固態堆芯多目標參數優化

4.1 固態堆芯優化問題

固態堆芯關鍵性能指標主要為堆芯最高溫度和堆芯內最大應力，其受到燃料棒、熱管之間排布間距以及各間隙大小等設計參數的影響。對固態堆芯尺寸進行參數優化，以期保持燃料總量不變，使優化后堆芯能獲得更低的堆內溫度和Mises應力。

基于ABAQUS、isight、機器學習等搭建聯合仿真智能優化平臺進行優化問題求解，優化技術路線如圖6所示。首先通過拉丁超立方抽樣在設計空間內構建樣本集，利用有限元求解堆芯溫度場和應力場，構成樣本數據庫，共獲得5 800組樣本數據。然后基于樣本數據庫，將5 000組用于機器學習訓練，800組用于測試，使用DNN建立固態堆芯設計參數與堆內最高溫度、最大Mises應力兩個關鍵性能指標的代理模型?；诖砟Ｐ?，使用NSGA-Ⅱ遺傳算法獲得連續的pareto解集，以獲取指定權重下的最優解。

圖6 優化流程技術路線圖Fig.6 Schematic diagram of optimization process

4.2 基于DNN構建代理模型

圖7為堆內最大Mises應力和最高溫度兩個優化目標的DNN預測值和有限元計算結果的對比。代理模型的擬合優度R2均大于0.998，且預測平均相對誤差均在0.8%以下，具有良好的預測準確性。

a——最大應力預測與有限元結果對比；b——最高溫度預測與有限元結果對比圖7 DNN代理模型預測與有限元計算結果對比Fig.7 Comparison of DNN surrogate model prediction and finite element results

4.3 多目標優化

利用NSGA-Ⅱ算法[27]對目標函數進行優化，圖8為堆內最大Mises應力和最高溫度的迭代曲線。經過5 000次迭代，搜集整個設計域中的非劣解集構成較為平滑的pareto前沿。圖9為最終得到的pareto解集。利用理想點法求解給定權重下的最優解，獲得了兩個優化目標權重相等時的最優解(240，967)。

圖8 應力、溫度迭代過程Fig.8 Iterative processes of stress and temperature

圖9 pareto解集與1∶1權重下的最優解Fig.9 pareto front solution and optimal solution with weight ratio of 1∶1

對等權重多目標優化的固態堆芯結構進行了有限元計算驗證，熱力耦合計算結果如圖10所示，其分布規律和初始設計十分相似。在表3中將優化后的堆芯結構和初始設計的各項參數進行對比，優化后堆內最高溫度為965 K，下降了8.44%，最大Mises應力為240 MPa，下降了34.43%。

圖10 初始設計與優化設計結果對比Fig.10 Result comparison between initial design and optimal design

表3 初始設計與優化設計參數對比Table 3 Parameter comparison between initial design and optimal design

5 結論

基于熱管冷卻反應堆固態堆芯工作特點，結合商用軟件ABAQUS和其子程序的二次開發，實現了對固態堆芯內含間隙的復雜傳熱行為的數值模擬以及腫脹、蠕變、熱應變等不同材料力學行為共同作用下堆芯力學行為數值模擬，得出結論如下。

1) 堆芯整體溫度分布具有顯著規律，燃料芯塊中心溫度較高，呈現沿著燃料芯塊徑向環狀下降的特點，此溫度場分布規律導致燃料芯塊中心受到向內的擠壓，而燃料芯塊外圍受拉，出現明顯的環狀拉應力區域，這些拉應力的產生是造成圓柱狀的燃料芯塊徑向開裂的重要原因。

2) 對固態堆芯尺寸設計參數進行優化，根據有限元大量計算獲得的數據樣本，使用DNN機器學習建立了尺寸參數和最高溫度、最大Mises應力的代理模型，使用NSGA-Ⅱ算法獲得了pareto前沿解集，以及在降低最高溫度和減小最大Mises應力兩個優化目標權重為1∶1的情況下最優解，相較于初始設計參數，最高溫度下降了8.44%，最大Mises應力下降了34.43%，改善了固態堆芯的性能。