韓梓超，章 靜,*，王明軍，田文喜，蘇光輝，秋穗正

(1.西安交通大學 能源與動力工程學院 動力工程多相流國家重點實驗室，陜西 西安 710049；2.西安交通大學 核科學與技術學院 陜西省先進核能技術重點實驗室，陜西 西安 710049)

核熱推進(NTP)是一種空間推進技術，其相關研究工作自20世紀40年代已經開始。核熱推進系統利用核反應堆釋放的能量直接加熱推進工質，高溫推進劑進入噴管中膨脹做功進而產生推力。在NASA的載人火星探測技術文件DRA5.0[1]中，核熱推進被列為首選的推進方案。核熱推進具有比沖高、推力大、壽命長等優(yōu)點[2]，是人類未來大規(guī)模空間活動的必然選擇[3]。研究表明，基于現有技術，在不大規(guī)模改動發(fā)動機尺寸的前提下，僅提高推進劑容量即可將往返火星的單程時間縮短一半[4]，可有效減少宇航員在宇宙中停留的時間，保障其安全。核熱推進是未來空間推進最有前景的動力形式，我國對于核熱火箭發(fā)動機的研制工作需要進一步開展[5]。

氫氣是空間核熱推進系統典型的推進劑，其流動換熱計算是核熱推進反應堆熱工水力計算的基礎。傳統的反應堆熱工水力計算多使用系統分析程序，如RELAP5、TRACE等。Wang等[6]使用RELAP5程序對CPR1000反應堆在全廠斷電事故下非能動余熱排出系統進行了分析計算，還對一個固有安全小型輕水堆在小破口事故下安全殼內各非能動系統的熱工水力特性進行了分析[7]，系統分析程序具有運算快、結果偏保守等優(yōu)點，在核反應堆研究領域得到了廣泛的應用。典型的適用于核熱推進系統的程序有ELM[8]和SAFSIM[9]等，ELM程序主要應用于核熱推進反應堆燃料元件穩(wěn)態(tài)熱工水力計算，Walton[8]使用該程序對NRX反應堆開展了相關計算，結果顯示流體溫度、壁面溫度、壓降等計算結果與實驗值符合良好，Walton在報告中還匯總了氫氣的流動換熱經驗關系式，具有較高的參考價值。國內有關核熱推進系統分析程序開發(fā)的工作較少，為切實保障核熱推進反應堆的研究設計工作，相關程序的開發(fā)應當盡早開展。

本研究基于現有的適用于壓水堆的系統分析程序進行二次開發(fā)，添加氫氣的物性模型及流動換熱和摩擦阻力關系式，使其適用于核熱推進系統內氫氣的流動換熱計算，采用公開文獻中的實驗數據及計算結果對相關模型進行驗證。此外，使用該二次開發(fā)的系統程序對一種典型的折流式核熱推進反應堆堆芯建模，并開展相應的計算及分析。

1 數值模型

1.1 氫氣物性模型

核熱推進系統中溫度、壓力變化范圍大，對氫氣物性模型的精度要求較高。本研究參考房玉良等[10]關于高溫氫介質物性的研究內容，針對核熱推進反應堆堆芯，所采用的氫氣物性模型適用范圍為溫度300～3 000 K、壓力0.1～15 MPa。比容、比熱容、熱導率和黏度的物性模型如下。

對于比容的求解采用Aungier-Redlich-Kwong(ARK)真實氣體模型[11]：

(1)

式中：a(T)、b、c均為系數，各系數的詳細表達式參見文獻[11]；p為壓力，Pa；R為氫氣氣體常數，R=4 124.5 J·kg-1·K-1；T為溫度，K；V為比容，m3·kg-1。

(2)

式中，ui、vi均為系數，各系數的值參見文獻[12]。實際氣體比定壓熱容與理想氣體比定壓熱容之間的關系式可參考文獻[10]。

對熱導率λ(ρ，T)的求解采用Assael模型[13]：

λ(ρ,T)=λo(T)+Δλ(ρ,T)+Δλc(ρ,T)

(3)

式中：λo(T)為稀薄氣體熱導率；Δλ(ρ，T)為附加熱導率；Δλc(ρ，T)為臨界增量；ρ為密度，kg·m-3。

對黏度η(T，ρ)的求解采用Muzny模型[14]：

(4)

式中：η0(T)為零密度極限條件下的黏度；η1(T)為相對于零密度條件由密度大于零引入的黏度增量；c1～c6為系數；Tr為對比溫度；ρr為對比密度。

1.2 氫氣流動換熱關系式

房玉良等[15]采用ANSYS FLUENT對高溫、高流速氫氣在圓管內的流動換熱特性做了研究，對比了程序計算值與經驗關系式求得的換熱系數，結果表明，修正的Miller-Taylor關系式符合較好，其表達式如下：

(5)

式中：Nu為努塞爾數；Re為雷諾數；Pr為普朗特數；L為通道長度；D為水力直徑；下標W、b表示壁面和主流。該關系式的適用范圍為：Re，3×104～4×105；TW/Tb，1.1～8。

1.3 氫氣摩擦阻力關系式

Walton在文獻[8]中整理了ELM程序采用的流動換熱關系式及摩擦阻力關系式，具有較高的參考價值。本研究模擬了氫氣在某燃料元件內的流動，對比了采用文獻[8]中不同摩擦阻力關系式進行計算的結果，具體內容見2.2節(jié)，在保證燃料元件壓降計算結果與參考值符合較好的前提下，推薦使用形式簡潔的Koo關系式：

(6)

式中，f為摩擦系數。

2 數值模型驗證

選用Taylor[16]的實驗數據驗證氫氣的物性模型，選用ELM程序對NRX反應堆燃料元件的模擬結果驗證流動換熱關系式和摩擦阻力關系式。

2.1 氫氣物性模型驗證

20世紀五六十年代，為研究核動力火箭發(fā)動機，美國開展ROVER項目，設計、建造并測試了多種氫氣冷卻反應堆及核動力火箭系統，期間開展了大量有關氫氣流動換熱的實驗。Taylor[16]發(fā)表了關于氫氣和氦氣在圓管內流動換熱實驗的研究成果，該實驗采用的試驗段由金屬鎢制成，試驗段長為0.227 m，內徑為2.946 mm，管壁厚為1.587 mm。該實驗測量了流量、進出口流體溫度和壓力以及管道外壁面溫度等參數，并整理得到摩擦系數和對流換熱系數的經驗關系式。

該報告中共有18組關于氫氣流動換熱的實驗數據，為充分驗證氫氣物性模型的正確性，本研究分別使用系統程序和ANSYS FLUENT對試驗段建模并計算，管道入口采用流量邊界，出口采用壓力邊界，給定入口流體的溫度和加熱功率。系統程序建模節(jié)點圖如圖1所示，入口處時間相關控制體給定入口流體溫度，時間相關接管給定入口流量，出口處時間相關控制體給定出口壓力，加熱功率由熱構件輸入。

圖1 圓管內氫氣流動換熱計算系統程序建模節(jié)點圖Fig.1 Nodalization of system code modeling of hydrogen flow and heat transfer calculation in circular tube

需要注意的是，文獻[16]中給出的溫度是總溫，而程序計算得到的溫度是靜溫，總溫與靜溫之間的關系為：

(7)

式中：Ttotal為總溫，K；Ts為靜溫，K；γ為比熱容比；Ma為馬赫數。

氫氣圓管內流動換熱計算結果與實驗值的對比如圖2所示，系統程序與CFD程序計算結果的相對誤差均在10%以內，而CFD程序的計算結果與實驗值符合更好。結果表明，采用1.1節(jié)所述的氫氣物性模型能較好地模擬氫氣在圓管內流動換熱的問題，出口總溫計算相對誤差小于10%，能滿足系統程序的計算要求，而對于簡單的氫氣流動換熱問題，CFD程序的模擬精度更高。

圖2 氫氣圓管內流動換熱計算結果與實驗值的對比Fig.2 Comparison of calculated result and experimental data of hydrogen flow in circular tube

2.2 流動換熱及摩擦阻力關系式驗證

Walton[8]使用ELM程序對NRX反應堆的六邊形燃料元件進行了計算，本研究選用NRX-A3反應堆燃料元件的實驗參數(表1)，燃料元件的幾何及簡化模型示意圖如圖3所示。

圖3 燃料元件幾何及簡化模型示意圖Fig.3 Scheme of geometry and simplified model of fuel element

表1 NRX-A3反應堆燃料元件實驗參數[8]Table 1 Experiment parameter of fuel element of NRX-A3 reactor[8]

使用系統程序對該燃料元件建模并計算，燃料元件簡化為單通道模型，模型、邊界條件與2.1節(jié)氫氣物性模型驗證部分相一致，節(jié)點圖如圖1所示，功率分布參考文獻[8]。采用不同摩擦阻力關系式進行計算，燃料元件通道內壓力的計算結果如圖4所示。取ELM程序計算結果的1%為相對誤差，結果顯示，采用Koo關系式、Taylor Ⅱ關系式和GRAPH-I-TITE G關系式計算得到的壓力結果均在誤差范圍內，而采用Taylor Ⅰ關系式和基本摩擦系數關系式計算得到的壓力在靠近燃料元件出口處超出了誤差范圍。綜合考慮結果的準確性和公式的簡潔性，推薦使用Koo關系式。

圖4 燃料元件壓力計算結果對比Fig.4 Comparison of pressure calculated result of fuel element

摩擦阻力關系式采用Koo關系式，流動換熱關系式采用修正的Miller-Taylor關系式，流體溫度及壁面溫度的計算結果與文獻值的對比如圖5所示。由圖5可看出，二次開發(fā)系統程序的計算結果與ELM程序的計算結果符合較好。除壓降、流體溫度和壁面溫度外，燃料元件內雷諾數與馬赫數的計算結果均符合良好，流動換熱及摩擦阻力關系式的選擇是合理的。

圖5 燃料元件溫度計算結果對比Fig.5 Comparison of temperature calculated result of fuel element

3 折流式核熱推進反應堆堆芯計算

3.1 折流式核熱推進系統簡介

折流即冷卻劑主流在堆芯不同區(qū)域內存在方向相反的流動，典型的折流式核熱推進系統如圖6所示，主要包括液氫儲箱、渦輪泵、再生冷卻通道、徑向反射層、外側環(huán)形堆芯、軸向反射層、內側柱形堆芯和噴管等部分，為方便表述，外側環(huán)形堆芯簡稱為外堆芯，內側柱形堆芯簡稱為內堆芯。低溫液氫在泵的作用下進入系統，隨后在噴管喉部再生冷卻通道中吸收熱量實現預熱，在此過程中，氫工質會跨越擬臨界點，物性變化比較劇烈[17]。超臨界氫氣隨后進入徑向反射層，吸收反射層釋放的熱量，之后推動渦輪，渦輪帶動泵轉動提供壓頭，此過程會伴隨較大的壓降。從渦輪流出的氫氣進入外堆芯，吸收外堆芯區(qū)域燃料元件釋放的熱量后，流經與外堆芯相鄰的軸向反射層，進入上腔室。在上腔室中，氫氣的流動方向發(fā)生改變，同時與來自不同燃料元件的氫氣充分交混，可有效降低氫氣溫度的不均勻程度[18]。交混后的氫氣流經與內堆芯相鄰的軸向反射層后，進入內堆芯進一步加熱至出口溫度達2 500 K以上。最后高溫氫氣在縮放噴管中膨脹做功，將內能轉化成動能，為飛行器提供推力。

圖6 折流式核熱推進系統示意圖Fig.6 Scheme of two-pass NTP system

3.2 反應堆堆芯幾何及建模

本研究針對的核熱推進反應堆采用CERMET燃料，將UO2或UN燃料顆粒彌散在鎢金屬基體中，最高運行溫度可達2 900～3 270 K，具有更安全、更經濟的優(yōu)勢[19]。若全堆采用金屬鎢作基體材料，由于密度較高，會導致堆芯質量過重，降低推重比。采用折流式反應堆設計，其堆芯結構如圖7所示，將核熱推進反應堆劃分為外堆芯和內堆芯兩個區(qū)域，溫度較低的外堆芯區(qū)域選用金屬鉬作基體材料，這樣可有效降低堆芯總質量?；w材料金屬鎢及金屬鉬的熱導率列于表2，密度及比定壓熱容列于表3。

表2 基體材料的熱導率[20]Table 2 Thermal conductivity of substrate material[20]

表3 基體材料的其他熱物理性質[20]Table 3 Other thermophysical property of substrate material[20]

圖7 折流式核熱推進反應堆堆芯示意圖Fig.7 Scheme of two-pass NTP reactor core

本文主要內容為氫氣相關模型的驗證及應用，運行參數、尺寸數據等僅根據經驗選取，初步設計堆芯高度為80 cm，總功率為495 MW，單個燃料元件的流通面積約為2.6×10-4m2。使用系統程序對該反應堆堆芯建模，為節(jié)省計算資源，取堆芯的1/6作為研究對象。程序采用單通道模型，將每個六邊形圈上的燃料元件等效為1根管，水力直徑取燃料元件上單個孔道的直徑，流通面積取所有孔道的流通面積之和，堆芯截面圖如圖8所示，其中不同顏色代表不同的單通道。

圖8 反應堆堆芯徑向截面示意圖Fig.8 Radial cross section scheme of nuclear reactor core

系統的入口采用流量邊界條件，同時給定入口流體的溫度；出口采用壓力邊界條件；堆芯燃料元件與軸向反射層相連，共等效為10根管型部件；上腔室等效為單一控制體；不同管型部件之間使用熱構件相連，模擬堆芯徑向不同位置的燃料元件之間的熱傳導，堆芯的釋熱功率通過熱構件輸入。系統建模的節(jié)點圖如圖9所示。

圖9 1/6堆芯建模節(jié)點圖Fig.9 Modeling nodalization of 1/6 reactor core

3.3 穩(wěn)態(tài)熱工水力計算

在前期的核熱推進反應堆堆芯設計工作中，通過合理設計不同燃料元件的富集度，可以展平徑向功率，不同通道的歸一化功率分布如圖10中柱狀圖所示，對于內堆芯區(qū)域的通道1～6，功率分布較為均勻，且沿徑向由內向外功率逐漸減小，對于外堆芯區(qū)域的通道7～10，通道7、8和9功率相近，通道10由于距離徑向反射層比較近，中子通量水平相對較高，功率較高。輸入某組設計工況的邊界條件進行計算，可得到相應的穩(wěn)態(tài)熱工水力結果，1/6堆芯的歸一化流量分布如圖10中散點圖所示。可以看出，由于內堆芯功率分布較為均勻，且氫氣在上腔室中混合后再流入內堆芯，內堆芯區(qū)域不同通道間的流量分布較為均勻。而外堆芯的流量分布不均勻程度較高，通道8和9由于功率水平較低，其流量高于外堆芯通道的平均值；特殊的是，對于與內堆芯區(qū)域相鄰的通道7，雖然其燃料元件功率低于外堆芯平均值，但其流量也低于平均值，且相差超過15%。

圖10 1/6堆芯歸一化功率及流量分布Fig.10 Normalized distribution of power and mass flow rate of 1/6 reactor core

該核熱推進反應堆采用的CERMET燃料為多孔棱柱形燃料，堆芯內的流動屬于典型的并聯通道間流動，并聯通道間的釋熱不均易引發(fā)流量分配的不均勻。由外堆芯區(qū)域功率及流量分布情況可看出，對于折流式核熱推進反應堆，功率無法真實反映釋熱的水平，需對比不同通道內流體的焓升來衡量堆芯釋熱的均勻程度。1/6堆芯歸一化焓升如圖11所示。通道6的焓升低于內堆芯的平均水平，而通道7的焓升高于外堆芯的平均水平，這是由于內外堆芯溫度差別較大，在交界處存在較強的熱傳導，部分內堆芯燃料元件釋放的熱量傳遞到外堆芯，因此在外堆芯區(qū)域，通道7的功率雖然較低，但其釋熱量遠高于平均值，在此影響下其流量遠低于平均值。

圖11 1/6堆芯歸一化焓升分布Fig.11 Normalized distribution of enthalpy rise of 1/6 reactor core

內外堆芯之間的熱傳導增加了堆芯釋熱的不均勻程度，造成個別燃料元件內的參數與平均值相差較大，危害系統的穩(wěn)定性。另外，對于外堆芯區(qū)域最內側的燃料元件(通道7)，其流量較低，釋熱量較高，易造成基體材料溫度超過限值，危害反應堆的安全。綜上所述，在折流式核熱推進反應堆堆芯的反應堆物理設計中，尤其對于外堆芯區(qū)域具有不同富集度的燃料元件的布置，由于內外堆芯之間存在較強的熱傳導，不能盲目地追求功率展平，應當結合堆芯的穩(wěn)態(tài)熱工水力計算，合理設計堆芯方案，使得冷卻劑流體溫度、流量等參數沿徑向分布盡量均勻，提高反應堆的安全性。

4 結論

本研究將氫氣的物性模型、流動換熱及摩擦阻力關系式植入一系統分析程序中，分別采用實驗數據及程序計算結果進行了驗證，使用該程序對折流式核熱推進反應堆堆芯建模并計算，所得結論如下。

1) 采用本文所述氫氣物性模型，對某圓管內氫氣流動換熱實驗的計算相對誤差小于10%；采用推薦的流動換熱及摩擦阻力經驗關系式，對某燃料元件內氫氣流動換熱的模擬結果與參考值符合較好。

2) 本研究采用的系統分析程序經二次開發(fā)，能應用于核熱推進系統的建模與計算，而程序的適用性、計算結果的準確性需要更多的驗證及分析。

3) 對于折流式核熱推進反應堆，內外堆芯之間的熱傳導增加了堆芯釋熱的不均勻程度，對堆芯穩(wěn)態(tài)特性影響較大，堆芯方案的設計應結合相應的熱工水力計算，確保溫度、流量等重要參數沿徑向分布均勻，提高安全性。