核熱推進(jìn)系統(tǒng)分析程序模型與計(jì)算方法初步研究

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

(清華大學(xué) 核能與新能源技術(shù)研究院,先進(jìn)反應(yīng)堆工程與安全教育部重點(diǎn)實(shí)驗(yàn)室,北京 100084)

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

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

1 NTP系統(tǒng)設(shè)計(jì)分析概念與內(nèi)涵

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

1.1 NTP系統(tǒng)設(shè)計(jì)分析思路

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

圖1 通用NTP設(shè)計(jì)程序框架Fig.1 Generic NTP design program framework

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

1.2 穩(wěn)態(tài)設(shè)計(jì)分析程序

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

除了NTP專用穩(wěn)態(tài)設(shè)計(jì)程序NESS外,NASA的GRC曾使用通用的NPSS(Numerical Propulsion System Simulation)程序?qū)TP系統(tǒng)進(jìn)行了穩(wěn)態(tài)性能分析工作。NPSS被航空航天工業(yè)廣泛用于葉輪機(jī)械、吸氣式推進(jìn)系統(tǒng)、液體火箭發(fā)動(dòng)機(jī)、發(fā)動(dòng)機(jī)控制系統(tǒng)和系統(tǒng)模型集成的建模[10]。該軟件具備初步設(shè)計(jì)、非設(shè)計(jì)點(diǎn)性能分析、瞬態(tài)性能分析和飛行試驗(yàn)數(shù)據(jù)關(guān)聯(lián)等能力[11]。

清華大學(xué)核能與新能源技術(shù)研究院自主開發(fā)了穩(wěn)態(tài)性能設(shè)計(jì)分析軟件PANES(Program for Analyzing Nuclear Engine Systems),該軟件具備NTP系統(tǒng)循環(huán)參數(shù)計(jì)算、顆粒床反應(yīng)堆輪廓方案、反射層冷卻方案、噴管構(gòu)型方案設(shè)計(jì)以及堆芯流動(dòng)換熱特性計(jì)算的能力[1,12]。

1.3 瞬態(tài)分析程序

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

表1 美國NTP瞬態(tài)分析程序Table 1 NTP transient analysis program of US

隨著計(jì)算技術(shù)發(fā)展與仿真技術(shù)的發(fā)展,瞬態(tài)分析軟件由專用簡化系統(tǒng)模型向通用分析平臺(tái)及全范圍系統(tǒng)模型過渡。而在建模方式上,這些系統(tǒng)分析程序依然采用了“流網(wǎng)-熱網(wǎng)”的建模思路。

2 NTP系統(tǒng)瞬態(tài)分析程序框架

按照NTP系統(tǒng)設(shè)計(jì)分析思路,為了形成一體化設(shè)計(jì)分析能力以提高設(shè)計(jì)效率,在現(xiàn)有穩(wěn)態(tài)程序PANES的基礎(chǔ)上采用“熱網(wǎng)-流網(wǎng)”的框架自主開發(fā)了NTP系統(tǒng)瞬態(tài)分析功能。

2.1 流網(wǎng)模型

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

(1)

(2)

(3)

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

圖2為流網(wǎng)基礎(chǔ)單元示意圖,選擇流網(wǎng)環(huán)節(jié)中進(jìn)口壓力、進(jìn)口比焓以及出口質(zhì)量流量為流網(wǎng)邊界,建立以“流阻/感-流容”(R/I-C)連接順序的流網(wǎng)基礎(chǔ)單元。多個(gè)基礎(chǔ)單元相連即可描述包含多個(gè)節(jié)點(diǎn)的任意流網(wǎng)環(huán)節(jié)。

圖2 流網(wǎng)基礎(chǔ)單元示意圖Fig.2 Schematic diagram of flow network basic unit

2.2 熱網(wǎng)模型

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

圖3 熱網(wǎng)基礎(chǔ)單元示意圖Fig.3 Schematic diagram of thermal network basic unit

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

(4)

(5)

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

2.3 瞬態(tài)求解算法

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

基于“流網(wǎng)-熱網(wǎng)”框架,可得到描述整個(gè)系統(tǒng)的數(shù)學(xué)模型,并建立系統(tǒng)性能分析的流程,如圖4所示。在穩(wěn)態(tài)分析中,可以將ODE中時(shí)間導(dǎo)數(shù)項(xiàng)置零,此時(shí),穩(wěn)態(tài)求解問題即轉(zhuǎn)化為非線性方程組求根問題,采用牛頓法即可解決。在瞬態(tài)分析中,由于隱式方法同樣需要求解非線性方程組,因此,在實(shí)現(xiàn)隱式方法的同時(shí)可以改進(jìn)穩(wěn)態(tài)分析程序,提高穩(wěn)態(tài)計(jì)算速度以實(shí)現(xiàn)快速的性能方案設(shè)計(jì)、分析以及優(yōu)化。高效的穩(wěn)態(tài)求解也可以為瞬態(tài)分析提供可靠的初始值。

圖4 NTP系統(tǒng)瞬態(tài)分析程序計(jì)算流程示意圖Fig.4 Schematic diagram of calculation process for transient analysis program in NTP system

3 反應(yīng)堆動(dòng)力學(xué)模型

反應(yīng)堆中子動(dòng)力學(xué)模型旨在解決NTP推力調(diào)節(jié)等瞬態(tài)過程中存在一系列反應(yīng)性引入以及反應(yīng)性反饋下的堆芯功率計(jì)算問題。在動(dòng)態(tài)問題中忽略堆芯內(nèi)中子密度隨空間分布可以引入點(diǎn)堆模型。

3.1 點(diǎn)堆模型動(dòng)態(tài)方程

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

(6)

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

其中:n為中子密度,cm-3;ρ為凈反應(yīng)性,$;Ci為先驅(qū)核濃度,cm-3;βi為第i組緩發(fā)中子份額;β為緩發(fā)中子總份額;λi為第i組衰變常數(shù),s-1;Λ為中子每代時(shí)間,s;I為緩發(fā)中子組數(shù)。

3.2 點(diǎn)堆模型求解流程

線性O(shè)DE的隱式求解需要在各時(shí)間步內(nèi)遞進(jìn)求解Ax=b,其計(jì)算流程如圖5所示。其中BE與CN方法的區(qū)別體現(xiàn)在如何構(gòu)造矩陣A和向量b上。在確定各時(shí)間步下凈反應(yīng)性ρ(t)后即可采用LU分解對線性方程組進(jìn)行求解。

圖5 點(diǎn)堆動(dòng)力學(xué)方程隱式求解流程圖Fig.5 Schematic diagram of implicit solution process for point reactor kinetics model

4 渦輪泵模型

渦輪泵是推進(jìn)劑供應(yīng)系統(tǒng)的核心部件,經(jīng)過噴管再生冷卻加熱后的推進(jìn)劑進(jìn)入渦輪做功驅(qū)動(dòng)離心泵工作以提供驅(qū)動(dòng)壓頭。渦輪泵瞬態(tài)分析模型是系統(tǒng)瞬態(tài)分析的基礎(chǔ),渦輪泵模型主要包括離心泵模型、渦輪模型以及轉(zhuǎn)軸模型3個(gè)部分。

4.1 離心泵模型

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

(7)

(8)

其中:Δp為揚(yáng)程,Pa;I為泵內(nèi)流體的慣性,kg/m4;ρR為流體的參考密度,kg/m3;LR為泵中流體的長度,m;AR為泵的參考流通面積,m2。I與泵的額定參數(shù)和泵實(shí)際的時(shí)間常數(shù)相關(guān)。

基于上述泵模型可獲得瞬態(tài)計(jì)算流程如圖6所示,采用該顯式計(jì)算方法需要知道當(dāng)前時(shí)刻參數(shù),以及下一時(shí)刻的邊界參數(shù)進(jìn)口參數(shù)、出口壓力、泵轉(zhuǎn)速,即可計(jì)算所需的泵瞬態(tài)參數(shù)。

圖6 離心泵瞬態(tài)計(jì)算流程Fig.6 Schematic diagram of transient calculation process for centrifugal pump

4.2 渦輪模型

(9)

(10)

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

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

4.3 轉(zhuǎn)軸模型

渦輪模型中假設(shè)離心泵葉輪、渦輪動(dòng)葉以及轉(zhuǎn)軸為剛體并繞轉(zhuǎn)軸中心定軸旋轉(zhuǎn)。離心泵和渦輪模型輸入?yún)?shù)均包括軸轉(zhuǎn)速,輸出參數(shù)均包括轉(zhuǎn)矩。通過式(11)剛體定軸轉(zhuǎn)動(dòng)角動(dòng)量方程即可將離心泵和渦輪進(jìn)行聯(lián)立瞬態(tài)求解。

(11)

其中:ω為渦輪泵軸轉(zhuǎn)速,rad/s;Jtp為渦輪泵轉(zhuǎn)動(dòng)部分轉(zhuǎn)動(dòng)慣量,kg·m2;Tt為渦輪驅(qū)動(dòng)轉(zhuǎn)矩,N·m;Tp為離心泵阻力力矩,N·m;Tf為摩擦阻力力矩,N·m。

5 結(jié)論

本文針對NTP系統(tǒng)研發(fā)需求,提出了系統(tǒng)性能設(shè)計(jì)主要內(nèi)容應(yīng)包括穩(wěn)態(tài)設(shè)計(jì)點(diǎn)性能分析與優(yōu)化、穩(wěn)態(tài)非設(shè)計(jì)點(diǎn)性能分析以及瞬態(tài)性能分析3個(gè)環(huán)節(jié)。為自主開發(fā)具備穩(wěn)態(tài)、瞬態(tài)設(shè)計(jì)分析能力的一體化的NTP系統(tǒng)分析程序,提出了基于“流網(wǎng)-熱網(wǎng)”的建模框架,建立了以常微分方程組為基礎(chǔ)的系統(tǒng)數(shù)學(xué)模型。基于上述框架,給出了點(diǎn)堆動(dòng)力學(xué)模型與渦輪泵瞬態(tài)模型的瞬態(tài)計(jì)算分析流程。后續(xù)將采用該框架進(jìn)一步完善系統(tǒng)與各部件模型,并對典型系統(tǒng)方案開展性能分析工作。