核熱推進系統分析程序模型與計算方法初步研究

毛晨瑞,吉 宇,孫 俊,郎明剛,石 磊

(清華大學 核能與新能源技術研究院,先進反應堆工程與安全教育部重點實驗室,北京 100084)

核熱推進(NTP)是利用核能將推進劑加熱到高溫從噴管膨脹并高速噴出產生推力的空間推進裝置,主要由推進劑供應系統、反應堆系統以及噴管等部件組成。推進劑供應系統主要包括推進劑貯箱、渦輪泵、管閥和測控儀表等部件,反應堆系統包括輻射屏蔽、堆芯和控制鼓等部件[1]。與化學火箭發動機相比,NTP具備更高的比沖和與之相當的推力。因此,在時效性要求較高的載人深空探測和軌道機動等任務方面具有巨大的應用潛力。美國和前蘇聯于20世紀50年代便開始圍繞洲際導彈推進、戰略導彈助推段攔截和大型有效載荷空間軌道轉移等場景開展了NTP系統研究,通過Rover/NERVA、Timberwind/SNTP和SEI等項目,取得了一系列重要成果,具體包括較為完備的NTP系統性能分析方法與程序、NTP反應堆方案、耐高溫抗氫蝕核燃料以及關鍵地面試驗設施等,前蘇聯更是建成了迄今為止成熟度最高的RD-0410核熱火箭發動機地面樣機[2-4]。但隨著化學推進的快速發展,NTP系統喪失了其原有的任務需求,眾多項目陸續中止。

近些年來,由于突出的性能優勢以及在眾多航天任務中的不可替代性,美國重新啟動了對NTP系統的研究,美國國家航空航天局(NASA)在2014年開始的核低溫推進級(NCPS)項目中以安全、高效、經濟作為主要的研發原則[5]。此外,美國國防先進研究計劃局(DARPA)在2020年牽頭實施的“敏捷的月空間飛行火箭驗證項目”(DRACO)目前正快速推進,并計劃于2027年實現在軌示范驗證[6]。由于NTP工作溫度極高,且反應堆系統與推進劑供應系統存在極強的耦合關系,存在實驗難度大和研發成本高等挑戰,所以針對NTP開發專用的系統設計分析程序是必要的。本文針對NTP系統性能分析需求,提出主要研發內容,并對已有的分析程序與框架進行評估,開展NTP系統瞬態模型以及計算方法研究。

1 NTP系統設計分析概念與內涵

NTP相比傳統的反應堆裝置存在以下差異:1) 系統為開放的一回路設計,規模小,部件間耦合緊密,動態響應快;2) 系統自身的設計裕量小,尤其是反應堆堆芯;3) 在快速啟動、停機和推力調節等典型瞬態過程中,存在系統和堆芯物理熱工參數的劇烈變化;4) 受空間應用的限制,無法設計類似于傳統反應堆的安全措施和冗余手段?；谝陨喜町?水堆或高溫氣冷堆所采用的系統設計與分析方法在NTP中可能會導致較大的偏差,甚至無法獲得合理可行的方案。因此有必要對NTP系統的設計與分析方法開展進一步研究。

1.1 NTP系統設計分析思路

根據NTP特點及前期探索,將系統設計和分析活動分為穩態設計點分析與優化、穩態非設計點性能分析與瞬態性能分析3部分,如圖1所示[7]。對NTP系統設計而言,首先需要選定裝置的一個主要技術狀態作為設計和分析的基準,也就是設計點,來進行熱力計算。在這個過程中,依據論證的系統構型進行熱力循環參數的確定和優化,具體涵蓋對要求的性能參數(主要是比沖和推力)、設計極限(如渦輪許用溫度限值和各部件的效率等)和設計選擇(如泵揚程、渦輪壓比和分流比等)之間的關系進行分析,獲得NTP系統性能參數的變化規律,作為優化的方向和依據。對空間裝置而言,有時還需要根據熱力循環參數對系統的質量進行初步估計,實現性能和質量的雙目標優化。在這個階段中,NTP系統可看作是柔性的“橡皮泥”,尺寸和參數都是不定的,有較大的自由度和優化空間。

圖1 通用NTP設計程序框架Fig.1 Generic NTP design program framework

經過設計點分析與優化確定初步的熱力循環方案后,則可以根據節點參數開展初步的部件和流道設計工作,獲得各部件的工作特性。非設計點性能分析是對處于NTP系統運行包線內的所有點進行分析,求解各部件的共同工作狀態,獲得系統性能參數的變化規律等。此時,可設想完成了一臺暫缺詳細結構信息的NTP裝置的“標準設計”。針對多種“標準設計”開展非設計點性能分析的工作,可進一步優化得到能滿足任務需求的NTP系統。非設計點性能分析后,可進一步研究系統和部件參數在各非設計點之間的動態變化過程,即瞬態性能分析。所獲得的特性可作為NTP系統啟動、停機和推力調節等過程運行方案設計的依據。同時瞬態性能分析工作也可用于研究反應性誤引入和葉輪機械斷軸等事故場景下的系統參數變化規律。對在動態過程中出現參數超過設計限值的情況,有必要與設計點優化和非設計點性能分析等工作迭代,從而給出更優的系統設計。

1.2 穩態設計分析程序

NTP系統設計點分析與優化以及非設計點性能分析均是基于穩態設計分析程序而開展的。NESS(Nuclear Engine System Simulation)是為了支持美國空間探索計劃(SEI)而由美國科學應用國際公司(SAIC)為NASA格倫研究中心(GRC)開發的NTP穩態設計分析程序。該程序結合了西屋公司開發的反應堆模型ENABLER與Aerojet通用公司開發的液體火箭發動機模擬程序ELES[8]。NESS可以設計和確定NTR推進系統組件(包括反應堆)的重量、性能和工作特性,以及發動機子系統參數。NESS程序與NERVA項目成果小型核火箭發動機(SNRE)設計結果進行了對比驗證[9]。

除了NTP專用穩態設計程序NESS外,NASA的GRC曾使用通用的NPSS(Numerical Propulsion System Simulation)程序對NTP系統進行了穩態性能分析工作。NPSS被航空航天工業廣泛用于葉輪機械、吸氣式推進系統、液體火箭發動機、發動機控制系統和系統模型集成的建模[10]。該軟件具備初步設計、非設計點性能分析、瞬態性能分析和飛行試驗數據關聯等能力[11]。

清華大學核能與新能源技術研究院自主開發了穩態性能設計分析軟件PANES(Program for Analyzing Nuclear Engine Systems),該軟件具備NTP系統循環參數計算、顆粒床反應堆輪廓方案、反射層冷卻方案、噴管構型方案設計以及堆芯流動換熱特性計算的能力[1,12]。

1.3 瞬態分析程序

表1列出了美國NTP瞬態分析程序情況,針對20世紀開展的NERVA、SNTP以及21世紀開展的NCPS等項目,美國開展了多型NTP系統瞬態分析程序的開發工作。國內的多家單位也開展了NTP系統設計、分析以及應用研究,劉忠恕[13]使用國產Modelica分析平臺MWorks建立了NTP系統整機模型,并開展了動態特性研究;鐘巖俊[14]使用AMESIM仿真程序開展了小推力核熱火箭的循環方案研究;Qi等[15]開發了NTPSS-vPower程序以開展NTP瞬態建模與研究。

表1 美國NTP瞬態分析程序Table 1 NTP transient analysis program of US

隨著計算技術發展與仿真技術的發展,瞬態分析軟件由專用簡化系統模型向通用分析平臺及全范圍系統模型過渡。而在建模方式上,這些系統分析程序依然采用了“流網-熱網”的建模思路。

2 NTP系統瞬態分析程序框架

按照NTP系統設計分析思路,為了形成一體化設計分析能力以提高設計效率,在現有穩態程序PANES的基礎上采用“熱網-流網”的框架自主開發了NTP系統瞬態分析功能。

2.1 流網模型

以“流阻/感-流容”的形式建立NTP系統分析的流網模型。噴管、閥門、離心泵、渦輪等部件模型可以簡化為單一的流阻/感模型。堆芯內流道以及再生冷卻通道需要在流阻/感模型的基礎上考慮流容模型。在流容模型中考慮流體的可壓縮性以及能量輸入,將質量連續性方程與能量守恒方程進行簡化后如式(1)、(2)所示。在流阻/感模型中將動量守恒方程簡化后可獲得考慮流體慣性以及阻力系數的質量流量計算方程,如式(3)所示。

(1)

(2)

(3)

其中:p為壓力,Pa;t為時間,s;ρ為流體密度,kg/m3;G為質量流量,kg/s;Kβ為體積彈性模量,Pa;a為流體聲速,m/s;V為控制體容積,m3;h為流體比焓,J/kg;Qi為流網節點與熱網節點間的換熱功率,W;Rf為阻力系數,Pa/(kg·s-1)2;A為通流面積,m2;L為管道長度,m;下標1、2分別表示對應控制體進、出口處參數。

圖2為流網基礎單元示意圖,選擇流網環節中進口壓力、進口比焓以及出口質量流量為流網邊界,建立以“流阻/感-流容”(R/I-C)連接順序的流網基礎單元。多個基礎單元相連即可描述包含多個節點的任意流網環節。

圖2 流網基礎單元示意圖Fig.2 Schematic diagram of flow network basic unit

2.2 熱網模型

與流網模型類似,以“熱容-熱阻”(R-C)的形式建立燃料組件、噴管壁等熱構件的熱網模型(圖3)。對于熱容環節,其連接的熱阻環節不局限于圖3中給出的形式,可根據建模的精度連接任意熱阻環節。式(4)為熱阻模型,用于描述固體-流體間對流換熱、固體與固體間的一維導熱,式(5)為能量守恒下描述熱構件溫度變化的熱容模型。

圖3 熱網基礎單元示意圖Fig.3 Schematic diagram of thermal network basic unit

q=A(T2-T1)/Rh=hA(T2-T1)

(4)

(5)

其中:q為熱流量,W;T為熱構件平均溫度,K;Rh為等效熱阻,(m2·K)/W;h為對流換熱系數,W/(m2·K);A為等效換熱面積,m2;m為熱構件質量,kg;c為熱構件平均比熱容,J/(kg·K)。

2.3 瞬態求解算法

通過建立“流網-熱網”的物理模型即可獲得描述NTP系統的數學模型即常微分方程組(ODE)。求解該方程組,便可得到NTP系統瞬態參數。對該問題,有顯式隱式兩類求解方式。與顯式方法相比,隱式方法通常更加穩定,具有更大的穩定時間步長范圍,對初始化要求更低。但是隱式方法中求解非線性方程組需要消耗更多計算資源。同時,由于隱式求解中,各部件關鍵參數需要耦合求解,提高了程序框架設計以及開發的難度。

基于“流網-熱網”框架,可得到描述整個系統的數學模型,并建立系統性能分析的流程,如圖4所示。在穩態分析中,可以將ODE中時間導數項置零,此時,穩態求解問題即轉化為非線性方程組求根問題,采用牛頓法即可解決。在瞬態分析中,由于隱式方法同樣需要求解非線性方程組,因此,在實現隱式方法的同時可以改進穩態分析程序,提高穩態計算速度以實現快速的性能方案設計、分析以及優化。高效的穩態求解也可以為瞬態分析提供可靠的初始值。

圖4 NTP系統瞬態分析程序計算流程示意圖Fig.4 Schematic diagram of calculation process for transient analysis program in NTP system

3 反應堆動力學模型

反應堆中子動力學模型旨在解決NTP推力調節等瞬態過程中存在一系列反應性引入以及反應性反饋下的堆芯功率計算問題。在動態問題中忽略堆芯內中子密度隨空間分布可以引入點堆模型。

3.1 點堆模型動態方程

在時空動力學的基礎上忽略空間效應即可獲得簡化的點堆模型,如式(6)所示,該模型的求解即為ODE的初值問題。由于點堆模型各組緩發中子的時間尺度差別較大,該常微分方程組具備較強的剛性。因此,采用后向歐拉(BE)以及Crank Nicholson(CN)等隱式求解方法可以采用更大的時間步長。

(6)

初值條件:n(t)|t=0=n0,Ci(t)|t=0=Ci0,i=1,2,…,I

其中:n為中子密度,cm-3;ρ為凈反應性,$;Ci為先驅核濃度,cm-3;βi為第i組緩發中子份額;β為緩發中子總份額;λi為第i組衰變常數,s-1;Λ為中子每代時間,s;I為緩發中子組數。

3.2 點堆模型求解流程

線性ODE的隱式求解需要在各時間步內遞進求解Ax=b,其計算流程如圖5所示。其中BE與CN方法的區別體現在如何構造矩陣A和向量b上。在確定各時間步下凈反應性ρ(t)后即可采用LU分解對線性方程組進行求解。

圖5 點堆動力學方程隱式求解流程圖Fig.5 Schematic diagram of implicit solution process for point reactor kinetics model

4 渦輪泵模型

渦輪泵是推進劑供應系統的核心部件,經過噴管再生冷卻加熱后的推進劑進入渦輪做功驅動離心泵工作以提供驅動壓頭。渦輪泵瞬態分析模型是系統瞬態分析的基礎,渦輪泵模型主要包括離心泵模型、渦輪模型以及轉軸模型3個部分。

4.1 離心泵模型

離心泵特性曲線由轉矩、壓頭關于流量系數曲線組成。泵的特性曲線實際上反映的是其穩態性能,在準穩態假設下可認為泵在不同工況下的性能是沿著其特性曲線變化的。式(7)給出了基于流網中“流阻/感”模型而建立的泵流量與增壓之間的關系:

(7)

(8)

其中:Δp為揚程,Pa;I為泵內流體的慣性,kg/m4;ρR為流體的參考密度,kg/m3;LR為泵中流體的長度,m;AR為泵的參考流通面積,m2。I與泵的額定參數和泵實際的時間常數相關。

基于上述泵模型可獲得瞬態計算流程如圖6所示,采用該顯式計算方法需要知道當前時刻參數,以及下一時刻的邊界參數進口參數、出口壓力、泵轉速,即可計算所需的泵瞬態參數。

圖6 離心泵瞬態計算流程Fig.6 Schematic diagram of transient calculation process for centrifugal pump

4.2 渦輪模型

(9)

(10)

氫渦輪模型具體的計算流程為:根據渦輪進口參數、出口背壓、轉速以及渦輪工作特性,迭代獲得渦輪內氣體質量流量和渦輪轉矩Tt,同時更新渦輪工作特性,如圖7所示。

圖7 渦輪模型計算流程圖Fig.7 Schematic diagram of turbine model calculation process

4.3 轉軸模型

渦輪模型中假設離心泵葉輪、渦輪動葉以及轉軸為剛體并繞轉軸中心定軸旋轉。離心泵和渦輪模型輸入參數均包括軸轉速,輸出參數均包括轉矩。通過式(11)剛體定軸轉動角動量方程即可將離心泵和渦輪進行聯立瞬態求解。

(11)

其中:ω為渦輪泵軸轉速,rad/s;Jtp為渦輪泵轉動部分轉動慣量,kg·m2;Tt為渦輪驅動轉矩,N·m;Tp為離心泵阻力力矩,N·m;Tf為摩擦阻力力矩,N·m。

5 結論

本文針對NTP系統研發需求,提出了系統性能設計主要內容應包括穩態設計點性能分析與優化、穩態非設計點性能分析以及瞬態性能分析3個環節。為自主開發具備穩態、瞬態設計分析能力的一體化的NTP系統分析程序,提出了基于“流網-熱網”的建模框架,建立了以常微分方程組為基礎的系統數學模型?；谏鲜隹蚣?給出了點堆動力學模型與渦輪泵瞬態模型的瞬態計算分析流程。后續將采用該框架進一步完善系統與各部件模型,并對典型系統方案開展性能分析工作。