鄧蛟龍，鄧　堅，李仲春，柴　翔,*

基于三維核熱耦合方法的吸氣式核反應堆瞬態(tài)特性研究

鄧蛟龍1，鄧堅2，李仲春2，柴翔1,*

（1. 上海交通大學，上海 200240；2. 中國核動力研究設計院，四川成都 610213）

吸氣式反應堆具有大比沖、高功率密度的優(yōu)良特點，在未來空天推進領域具有廣闊的應用前景。該反應堆在正常運行條件下需要完成各類功率快速響應的瞬態(tài)工況，并在極短時間內(nèi)完成啟停堆過程，采用傳統(tǒng)的數(shù)值方法研究其瞬態(tài)三維核熱耦合特性存在較大的困難。本文基于開源OpenFOAM平臺，開發(fā)了三維瞬態(tài)核熱耦合求解程序，提出了新型流固熱耦合邊界，實現(xiàn)了固體堆芯與流道內(nèi)部流動傳熱特性的快速預測，使用點堆模型求解堆芯中子物理特性，通過耦合物理模塊與熱工模塊實現(xiàn)了瞬態(tài)過程中氣式反應堆堆芯流動傳熱特性的數(shù)值分析。模擬結果表明，吸氣式反應堆在功率快速提升的同時，固體材料的溫度不會迅速升高，而是平穩(wěn)上升直至穩(wěn)態(tài)，反應堆在瞬時熱功率超過設計基準值的啟堆過程中仍然是安全的。

流固耦合；核熱耦合；OpenFOAM

隨著人類空間探索進程的不斷深入，化學推進技術已經(jīng)逐漸成為制約航天器飛行壽命與航行距離的重要因素，吸氣式反應堆為突破這種制約提供了可能性，在空天推進領域具有較高的應用價值。美蘇在20世紀開展過多次針對吸氣式反應堆發(fā)動機的實驗研究，主要包括核渦輪發(fā)動機XNJ-140E[1]以及核沖壓發(fā)動機TORY系列[2]等。然而該反應堆具有技術難度大、研發(fā)周期長的特點，兩國歷經(jīng)多年研究仍未達到實際應用的階段。我國在這一領域的研究尚處于理論驗證的起步階段，對于該反應堆安全方面的設計可行性仍存在疑問，需要進行大量的數(shù)值模擬研究為后續(xù)實驗提供支撐。該發(fā)動機結構如圖1所示。

圖1　吸氣式發(fā)動機整體結構

空天推進的特殊應用背景要求吸氣式反應堆工作在高溫、高壓、高功率密度的極限工況，以盡可能增大比沖，提高推進效率。相較于傳統(tǒng)反應堆數(shù)小時的啟堆時間，吸氣式反應堆的一大特點在于通過精準的反應性控制機制實現(xiàn)的功率快速變化，啟堆往往在數(shù)十秒內(nèi)完成，需要評估該過程中堆芯中子與熱工參數(shù)的反饋效應。因此針對這一反應堆的數(shù)值模擬是強瞬態(tài)工況下的核熱耦合過程。

本文針對吸氣式反應堆多流道、固體燃料結構與高溫空氣冷卻的特點，探究核熱耦合方法并利用開源CFD軟件OpenFOAM[3]進行耦合平臺的開發(fā)，利用該平臺實現(xiàn)對吸氣式反應堆的三維全堆芯瞬態(tài)響應數(shù)值模擬。文中工作將對該反應堆的后續(xù)安全設計提供一定的參考。

1　數(shù)值模型

1.1　模擬系統(tǒng)

以TORY反應堆的設計方案為基礎對吸氣式反應堆進行優(yōu)化建模，優(yōu)化后的反應堆堆芯關鍵參數(shù)如圖2與表1所示。冷卻劑流道位于六邊形固體燃料元件中心位置，整個堆芯由燃料元件緊密排布構成，堆芯軸向與徑向外部被反射層包圍，上下反射層具有不同的厚度。

反應堆冷卻劑為高溫空氣，經(jīng)堆芯底部流入頂部流出。燃料是以BeO為基體的93.5%富集度UO2陶瓷基體燃料，鈾元素均勻彌散在BeO基體中。此外，燃料中混合有少量的ZrO2與Y2O3，其中ZrO2用于減輕高溫工質(zhì)對固體結構的腐蝕；Y2O3用于減輕高溫條件下鈾氧化物的揮發(fā)。反射層材料為BeO，用于減輕堆芯的中子泄露。模擬中各區(qū)域的相關物性來源于文獻[4]。

結構上反應堆具有對稱性，為節(jié)省計算時間，模擬區(qū)域選定為圖2（b）所示的不含上下反射層區(qū)域的十二分之一堆芯。模擬區(qū)域具有大量冷卻劑流道結構，孔隙率達到53%。使用傳統(tǒng)的CFD方法研究流固結構之間的瞬態(tài)傳熱問題存在較大困難。

表1　堆芯關鍵參數(shù)

圖2　吸氣式反應堆堆芯結構

1.2　中子學模型

目前針對中子輸運方程的求解主要采用確定論方法或蒙特卡羅方法。依據(jù)輸運方程在空間、角度和能量上不同程度的簡化，確定論方法進一步分為點堆動力學模型、中子擴散模型和中子輸運模型[5]。本文選用點堆模型計算反應性與功率之間隨時間變化的關系，公式如下：

總反應性可表示為下式：

由于氣體冷卻介質(zhì)的密度較低，因此冷卻劑密度效應暫不在本文中考慮。

為了獲取吸氣式反應堆關鍵中子學參數(shù)，本文采用蒙特卡羅程序OpenMC[7]對該反應堆進行建模以得到不同燃料溫度下反應堆的多普勒反應性，擬合得多普勒反應性的計算公式如下：

在功率提升的瞬態(tài)過程中，外部反應性通過人為給定一線性關系，如公式（5）：

圖3　外部反應性變化

Fig.3　The change of external reactivity

點堆模型無法描述中子通量的空間分布特性。為了在計算中考慮功率的空間分布，本文基于OpenMC的計算結果重構出堆芯功率的空間分布特性，如圖4所示。

圖4　歸一化功率分布

1.3　堆芯流動傳熱模型

在這種計算方法中，流固交界兩側網(wǎng)格越精細，溫度梯度的計算結果越精確，更能準確預測壁面處的熱流密度。為了精確模擬流體區(qū)域近壁面的流動傳熱狀況，通常在近壁面處劃分大量流體網(wǎng)格，由于固體堆芯中存在近千根冷卻劑通道，使用這種傳統(tǒng)流固耦合模型在有限的計算資源下無法對吸氣式反應堆進行全堆芯模擬計算。

本文提出一種改進的流固熱耦合邊界模型用于減少流體區(qū)域的網(wǎng)格數(shù)量，使用牛頓冷卻定理計算流體側熱流密度，即將公式（11）更改為下式：

2　吸氣式反應堆瞬態(tài)工況研究

2.1　工況設置

基于本文開發(fā)的流固耦合傳熱計算模型與多孔介質(zhì)模型對吸氣式核反應堆數(shù)值模擬的簡化，數(shù)量龐大的冷卻劑流道無需逐一建立模型與劃分網(wǎng)格。冷卻劑區(qū)域作為整體進行模擬，極大減少冷卻劑區(qū)域的網(wǎng)格數(shù)量。流體幾何與固體燃料幾何基本重合，燃料幾何將保留各圓形冷卻劑流道孔洞。為合理調(diào)配計算資源，需對網(wǎng)格數(shù)量進行敏感性分析，設置的4套網(wǎng)格方案如表2所示。燃料區(qū)域網(wǎng)格數(shù)量最多，對計算精度的影響相對較大，因此前3組網(wǎng)格用于分析燃料區(qū)域網(wǎng)格數(shù)量；第1組和第4組網(wǎng)格用于分析冷卻劑區(qū)域網(wǎng)格數(shù)量的影響。以表3所示的滿功率穩(wěn)態(tài)工況作為邊界條件進行網(wǎng)格敏感性分析，計算獲得各網(wǎng)格穩(wěn)態(tài)時堆芯燃料軸向壁面溫度隨高度分布如圖5（a）所示。綜合考慮計算精度與計算資源消耗之后，本文選取網(wǎng)格1進行后續(xù)數(shù)值模擬，其二維截面結構如圖5（b）所示。

表2　網(wǎng)格數(shù)量

圖5　網(wǎng)格敏感性分析

2.2　模擬結果

穩(wěn)態(tài)滿功率計算結果如圖6所示，圖6（a）為功率分布，本文模擬中基于圖4所示的功率分布來考慮熱功率的空間分布特性。圖6（b）為燃料的溫度分布，圖6（c）、（d）分別為冷卻劑的速度和溫度分布。由于空間分布熱源的影響，溫度速度等熱工水力參數(shù)呈現(xiàn)相似的分布特征。冷卻劑自下向上流動，因此堆芯上半部分區(qū)域換熱效果較差，各熱工參數(shù)總體上在堆芯靠近出口處的中心位置出現(xiàn)較高值。

為了定量分析滿功率穩(wěn)態(tài)時燃料區(qū)域的溫度熱點位置，現(xiàn)取圖6（b）中兩條特征線上的冷卻劑與燃料的溫度數(shù)據(jù)分別繪制如圖7（a）、（b）所示。穩(wěn)態(tài)時燃料與冷卻劑區(qū)域已經(jīng)建立起較為穩(wěn)定的溫差，從能量守恒的角度，徑向溫差不斷減小，軸向溫差先增大再減小，趨勢與能量空間分布一致。從圖7（c）有，沿軸向上燃料溫度在不斷增大，靠近出口處的燃料溫度增速減緩，并在堆芯近出口的徑向中心處取到全堆芯燃料溫度的最大值1 604 K。

表3　模擬工況設置

圖6　工況1三維分布

圖6　工況1三維分布（續(xù)）

圖7　工況1定量分析

圖8　工況4功率提升2倍計算結果

圖9　工況5～8啟堆過程計算結果

3　結論

本文利用OpenFOAM針對吸氣式反應堆開展三維全堆芯瞬態(tài)核熱耦合模擬，獲得的主要結論如下：

（1）開發(fā)了改進的流固熱耦合邊界模型，求解得到的流體主流溫度滿足能量守恒方程，求得的固體壁面溫度取決于傳熱關系式的精度。該方法可以有效降低流體區(qū)域的網(wǎng)格數(shù)量。

（2）開發(fā)了瞬態(tài)核熱耦合求解器，以燃料的多普勒效應作為耦合的關鍵接口。該求解器能夠?qū)崿F(xiàn)同一套網(wǎng)格體系下的內(nèi)耦合瞬態(tài)求解。

（3）對反應堆的滿功率穩(wěn)態(tài)模擬表明，燃料溫度熱點出現(xiàn)在靠近出口的堆芯中心位置，約為1 604 K，不會發(fā)生水侵事故。

（4）吸氣式反應堆瞬態(tài)啟堆工況模擬表明，功率提升呈現(xiàn)階段性的劇烈變化，應該提供更精確的反應性控制機制緩解功率在反應性轉(zhuǎn)折點處發(fā)生的突變現(xiàn)象。受到堆芯熱容的影響，燃料溫度卻不會迅速升高，將平滑上升直至穩(wěn)態(tài)。這一發(fā)現(xiàn)對該反應堆的安全設計具有重要參考價值。

［1］ Chase C L.APEX-908 XNJ140E Nuclear Turbojet Part B［R］．Flight Propulsion Laboratory Department，1962：17-22.

［2］ Hadley J W.Tory II-A a Nuclear Ramjet Test Re-actor［R］．U.S.Atomic Energy Commission，1959：20-21.

［3］ Moukalled F，Mangani L，Darwish M.The Finite Volume Method in Computational Fluid Dynamics［M］．Switzerland：Springer International Publishing，2016.

［4］ Amaral J.Thermophysical Properties of Materials for Nuclear Engineering a Tutorial and Collection of Data［M］．Vienna：IAEA，2009.

［5］ Wang J，Wang Q，Ding M.Review on Neutronic/Thermal- hydraulic Coupling Simulation Methods for Nuclear Reactor Analysis［J］．Annals of Nuclear Energy，2019，137.

［6］ Sun H. et al. Conceptual Design and Analysis of aMultipurpose Micro Nuclear Reactor Power Source［J］．Annals of Nuclear Energy，2018，121：118-127.

［7］ Yu J. Chapter 6-Nuclear Power Plant Design an-d Analysis Codes［M］．Sawston：Woodhead Publishing，2021：123-138.

［8］ Tuominen R. Coupling Serpent and OpenFOAM f-or Neutronics-CFD Multiphysics Calculations［D］．A-alto University，2015：22-23.

［9］ Hilsenrath，J. C. Tables of Thermal Properties of G-ases［M］．NBS Circular，1955，564.

［10］ Taylor M F. Prediction of Friction and Heat Tra-nsfer Coefficients with Large Variations in Fluid Properties［R］．NASA Lewis Research Center，1970：5-7.

［11］ Walter C E. Engineering Design of The Tory II-C Nuclear Ramjet Reactor［R］．University of Calif-ornia，1964：2-10.

Study on Transient Characteristics of Nuclear Ramjet Reactor Based on Three-Dimensional Nuclear Thermal Coupling Method

DENG Jiaolong1，DENG Jian2，LI Zhongchun2，CHAI Xiang1,*

（1. Shanghai Jiao Tong University，Shanghai 200240，China；2. Nuclear Power Institute of China，Chengdu of Sichuan Prov. 610213，China）

The ramjet reactor has excellent characteristics of large specific impulse and high power density，and has broad application prospects in the field of aerospace propulsion in the future. Under normal operating conditions，the reactor needs to complete various transient operating conditions with rapid power response，and complete the startup and shutdown process in a very short time. The use of traditional numerical methods to study its transient three-dimensional nuclear-thermal coupling characteristics has great difficult. Based on the open source OpenFOAM platform，this paper developed a three-dimensional transient nuclear-thermal coupling solution program，proposed a new fluid-solid coupling boundary，and realized the rapid prediction of the flow and heat transfer characteristics of the solid core and the flow channel. The point kinetic model was employed to solve the core neutron physical properties.Through the coupling of the physical module and the thermal module，realize the numerical analysis of the flow and heat transfer characteristics of the gas reactor core during the transient process. The simulation results show that while the power of the ramjet reactor increases rapidly，the temperature of the solid material will not rise fast，but will rise steadily until it reaches a steady state. The reactor is still safe during the startup process when the instantaneous thermal power exceeds the design benchmark value.

Fluid-solid coupling；Nuclear-thermal coupling；OpenFOAM

TL48

A

0258-0918（2022）01-0018-10

2021-07-06

國家自然科學基金（51806139，11922505）；上海市工業(yè)強基專項資助（GYQJ-2018-2-02）

鄧蛟龍（1997—），男，四川攀枝花人，在讀研究生，現(xiàn)主要從事熱工水力方面研究

柴翔，E-mail：xiangchai@sjtu.edu.cn