脈沖反應(yīng)堆脈沖運(yùn)行瞬態(tài)的子通道分析

陳立新，朱 磊，馬騰躍，江新標(biāo)，陳 達(dá)

（1.南京航空航天大學(xué) 核科學(xué)與工程系，江蘇 南京 210016；2.西北核技術(shù)研究所，陜西 西安 710024）

脈沖反應(yīng)堆是一種以鈾氫鋯為燃料的池式研究堆，由于鈾氫鋯材料具有較大的瞬發(fā)負(fù)溫度系數(shù)，因此該類型反應(yīng)堆既可穩(wěn)態(tài)運(yùn)行，也可脈沖運(yùn)行。對(duì)于脈沖運(yùn)行工況，反應(yīng)堆的功率變化快，脈沖功率峰的寬度達(dá)ms量級(jí)，功率變化的范圍大，脈沖前的穩(wěn)態(tài)功率與脈沖峰值功率相差6～7個(gè)量級(jí)，這種快速的功率變化對(duì)堆芯的傳熱安全提出了更高要求。準(zhǔn)確計(jì)算脈沖運(yùn)行狀態(tài)下堆芯的熱工參數(shù)是評(píng)價(jià)脈沖反應(yīng)堆瞬態(tài)安全的重要內(nèi)容。

子通道模型是目前核反應(yīng)堆堆芯熱工水力分析中較精細(xì)的計(jì)算模型。由于子通道模型能準(zhǔn)確給出反應(yīng)堆各組件或各燃料通道內(nèi)的流量、溫度和壓力等熱工參數(shù)的分布，故基于子通道模型的分析程序成為反應(yīng)堆熱工水力分析的重要工具。世界上主要核能國家分別推出了各系列的子通道分析程序，如THINC 系列［1］、COBRA 系列［2］等。這些較成熟的商用子通道程序大多針對(duì)動(dòng)力反應(yīng)堆設(shè)計(jì)，其適用范圍也多集中在高溫、高壓、高質(zhì)量流量的強(qiáng)迫循環(huán)冷卻方式，對(duì)脈沖堆這類低溫、常壓、低質(zhì)量流量的自然循環(huán)冷卻反應(yīng)堆缺乏適用性。由于以往的脈沖堆瞬態(tài)熱工安全分析多采用較保守的單通道方法，無法準(zhǔn)確給出堆芯內(nèi)部詳細(xì)的熱工安全參數(shù)的分布，因此，本文基于子通道分析方法，開發(fā)脈沖堆瞬態(tài)特性分析的子通道程序，并以西安脈沖反應(yīng)堆為研究對(duì)象，分析堆芯的瞬態(tài)熱工安全特性。

1 子通道物理模型

1.1 子通道模型的數(shù)學(xué)方程

一般，子通道分析常采用均勻平衡模型，即假設(shè)氣、液兩相的相對(duì)速度為零，且具有相同的壓力與溫度?；谝陨霞僭O(shè)，可將冷卻劑看作氣、液兩相的混合物而進(jìn)行統(tǒng)一處理，其連續(xù)性方程及能量和動(dòng)量方程［3］如下。

連續(xù)性方程為：

能量守恒方程為：

動(dòng)量守恒方程為：式中：ρ為流體密度；τ 為時(shí)間；u 為速度矢量；ui、uk（i、k=1，2，3）分別為不同的速度分量；E 為流體具有的能量；T 為應(yīng)力張量；f 為體積力；Q 為流體發(fā)熱量。

以上方程組的求解非常困難，為簡化方程，一般的子通道分析方法均需進(jìn)行必要假設(shè)，認(rèn)為通道中冷卻劑的軸向流速度遠(yuǎn)大于橫流速度，以至于橫流在通過通道間的間隙后就失去了其原來具有的方向而隨軸流運(yùn)動(dòng)。這樣就可將動(dòng)量方程分解為軸向動(dòng)量守恒方程和橫向動(dòng)量守恒方程。以上子通道方程的推導(dǎo)參見文獻(xiàn)［3］。

1.2 子通道方程的數(shù)值解法

對(duì)于子通道方程，可采用有限差分方法對(duì)其進(jìn)行求解。將堆芯劃分為N 個(gè)子通道，并將各通道在軸向上分為L 層，將每層看作1個(gè)控制體，如圖1所示。在控制體內(nèi)部，流體具有單一的溫度、壓力、密度、比容等參數(shù)。將子通道守恒方程應(yīng)用于每個(gè)控制體，則可得到N×L組守恒方程。

對(duì)于離散化方程守恒方程，在每個(gè)時(shí)間步長Δτ內(nèi)按軸向分層逐層對(duì)各通道進(jìn)行迭代求解，得到本層的所有參數(shù)后，進(jìn)入下一層的求解，其中，能量方程需同燃料元件導(dǎo)熱方程聯(lián)合求解。

圖1 子通道網(wǎng)格劃分Fig.1 Mesh generation of sub-channel

2 輔助模型

2.1 功率計(jì)算模型

瞬態(tài)分析時(shí)，由于功率是時(shí)間的函數(shù)，且堆芯溫度的改變會(huì)影響功率的變化，因此需耦合熱工分析模塊與功率計(jì)算模塊進(jìn)行計(jì)算。為避免復(fù)雜的物理熱工耦合，本文采用考慮溫度反饋的點(diǎn)堆動(dòng)態(tài)方程計(jì)算堆芯的功率變化，同時(shí)也可根據(jù)物理計(jì)算程序的功率模擬結(jié)果將各燃料元件或堆芯整體的功率-時(shí)間變化曲線作為外部參數(shù)輸入。

2.2 燃料元件導(dǎo)熱模型

對(duì)于燃料元件的導(dǎo)熱，采用粗棒非穩(wěn)態(tài)一維導(dǎo)熱模型，該模型可用于模擬燃料元件穩(wěn)態(tài)和瞬態(tài)導(dǎo)熱問題，其通用控制方程為：

式中：r為與熱量傳遞方向平行的坐標(biāo)；F（r）為與導(dǎo)熱面積有關(guān)的因子；Q 為內(nèi)熱源；λ為導(dǎo)熱系數(shù)；T 為溫度。對(duì)于式（4），可將其在時(shí)間和空間上進(jìn)行離散，采用差分法進(jìn)行求解。

2.3 熱傳輸模型

對(duì)于燃料元件包殼與冷卻劑之間的對(duì)流換熱系數(shù)的計(jì)算，采用對(duì)流換熱的分區(qū)處理方法，即傳熱曲線被劃分為4個(gè)分區(qū)［4］。單相對(duì)流傳熱區(qū)采用Sieder-tate關(guān)系式，過冷沸騰區(qū)采用McAdams關(guān)系式，過渡沸騰傳熱采用對(duì)數(shù)坐標(biāo)的沸騰曲線，取臨界熱流密度和最小熱流密度線性內(nèi)插法估算傳熱系數(shù)，膜態(tài)沸騰傳熱采用簡化的處理方法。圖2為本文計(jì)算的脈沖堆傳熱曲線。

圖2 脈沖堆燃料包殼與冷卻劑的傳熱曲線Fig.2 Heat transfer curve of pulsed reactor between clad and coolant

2.4 水和水蒸氣物性的計(jì)算模型

選用水和水蒸氣性質(zhì)國際聯(lián)合會(huì)（IAPWS）發(fā)布的IAPWS-IF97公式［5］計(jì)算水和水蒸氣的物性。IAPWS-IF97公式將水和蒸氣按性質(zhì)的不同分成5個(gè)區(qū)：1區(qū)為常規(guī)水區(qū)，2區(qū)為過熱蒸氣區(qū)，3區(qū)為臨界區(qū)，4區(qū)為濕蒸氣區(qū)，5區(qū)為低壓高溫區(qū)。本文主要用到1、2 區(qū)的數(shù)據(jù)。IAPWS-IF97 適 用 范 圍 為273.15 K ≤T ≤2 273.15K，壓力p≤100MPa。每個(gè)分區(qū)的物性參數(shù)計(jì)算公式參見文獻(xiàn)［5］。

3 計(jì)算程序的編制及驗(yàn)證

3.1 程序編制

在以上對(duì)子通道模型進(jìn)行離散化的基礎(chǔ)上，結(jié)合燃料元件導(dǎo)熱模型、包殼冷卻劑對(duì)流換熱模型、瞬態(tài)功率計(jì)算模型、水及水蒸氣物性計(jì)算模型等編制脈沖堆瞬態(tài)安全分析的子通道程序PRC-STAC。程序采用FORTRAN 語言編寫，內(nèi)部計(jì)算模塊采用子函數(shù)形式，主程序通過輸入輸出控制對(duì)各子函數(shù)的調(diào)用。程序主要包括輸入輸出模塊、功率計(jì)算與處理模塊、燃料元件溫度場分析模塊、換熱系數(shù)計(jì)算模塊、水力計(jì)算模塊、輔助參數(shù)計(jì)算模塊等功能模塊。

3.2 程序驗(yàn)證

本文通過對(duì)OSTR 堆［6］典型的瞬態(tài)計(jì)算結(jié)果開展比對(duì)計(jì)算，來驗(yàn)證PRC-STAC程序的可靠性。OSTR 反應(yīng)堆屬于美國GA 公司研制的TRIGA MARK Ⅱ型反應(yīng)堆，與脈沖反應(yīng)堆結(jié)構(gòu)類似。該堆額定運(yùn)行功率為1.1 MW，反應(yīng)堆工作于常壓環(huán)境，堆芯依靠水的自然循環(huán)冷卻。本文選取了RELAP 5-3D 軟件計(jì)算的OSTR 反應(yīng) 堆 采 用LEU 燃 料（235U 富 集 度 為19.75%，與脈沖堆相同）的堆芯瞬態(tài)計(jì)算結(jié)果。為使二者更具可比性，選擇文獻(xiàn)［6］中8通道計(jì)算模型給出的計(jì)算結(jié)果進(jìn)行比對(duì)。

圖3 示出了兩個(gè)程序計(jì)算的OSTR 堆堆芯燃料最高溫度隨時(shí)間的變化曲線計(jì)算比對(duì)結(jié)果。計(jì)算的燃料溫度變化趨勢(shì)一致，燃料到達(dá)最高溫度的時(shí)間略有差別。圖4為燃料溫度達(dá)最大值時(shí)溫度在燃料棒徑向上的分布計(jì)算結(jié)果。由比對(duì)結(jié)果可看出，本文的計(jì)算結(jié)果與RELAP5-3D的符合較好。

圖3 燃料最高溫度隨時(shí)間的變化Fig.3 Fuel maximum temperature vs time

4 脈沖堆瞬態(tài)特性分析

為準(zhǔn)確模擬反應(yīng)堆在脈沖工況下堆芯熱工參數(shù)的時(shí)間-空間變化，采用三維時(shí)空動(dòng)力學(xué)程序XAPR-HENKO 對(duì)脈沖堆在引入3.5＄脈沖時(shí)的堆芯三維功率分布進(jìn)行計(jì)算［7］，脈沖工況的動(dòng)力學(xué)參數(shù)列于表1。

圖4 燃料元件溫度沿徑向的變化Fig.4 Fuel element temperature variation along radial direction

表1 脈沖工況的動(dòng)力學(xué)參數(shù)Table 1 Kinetics parameters in pulse operation

分析時(shí)選取沿堆芯徑向布置的1組燃料元件進(jìn)行計(jì)算，其柵元編號(hào)為D2、E2、F2、G2、H2、I2，其在脈沖堆芯內(nèi)的位置如圖5 所示。圖6為XAPR-HENKO 程序計(jì)算的以上幾根燃料元件的功率隨時(shí)間的變化曲線。

4.1 單棒計(jì)算

D2測(cè)溫燃料元件位于脈沖堆芯脈沖運(yùn)行時(shí)的熱點(diǎn)位置，因此首先對(duì)D2元件及其周圍的冷卻劑進(jìn)行分析。圖7 為發(fā)射3.5＄脈沖后不同時(shí)刻D2 燃料元件內(nèi)部燃料溫度在不同時(shí)刻的徑向分布。從圖7可看出，在發(fā)射脈沖的最初階段，燃料溫度的分布呈邊緣高、中間低的分布特征，這主要是由于粗棒元件自屏效應(yīng)的影響。隨時(shí)間的變化，燃料溫度逐漸趨于一致。

圖5 選取的燃料元件在脈沖堆芯的位置Fig.5 Position of fuel elements in core

圖6 燃料元件功率隨時(shí)間的變化Fig.6 Fuel element power vs time

圖7 脈沖時(shí)D2元件燃料徑向溫度變化Fig.7 Fuel radial temperature variation of D2element in pulse operation

由于在脈沖發(fā)射時(shí)，燃料邊緣的溫度最高，圖8示出了燃料邊緣和包殼溫度隨時(shí)間的變化曲線。

脈沖進(jìn)行時(shí)，燃料邊緣的溫度隨反應(yīng)堆功率的上升迅速升高，而包殼溫度的上升則滯后于燃料溫度的變化，顯然，這是由于熱量在向包殼傳遞過程中存在時(shí)間延遲。燃料內(nèi)部的能量釋放來自核裂變反應(yīng)，可認(rèn)為瞬時(shí)釋放全部能量，熱量在短時(shí)間內(nèi)來不及向包殼傳遞；因此包殼的溫度變化會(huì)滯后于燃料邊緣，同時(shí)由于包殼與冷卻劑直接接觸，其溫度的最大值也較燃料溫度低很多，但最高溫度也接近500 ℃。這一瞬間高溫持續(xù)時(shí)間極短，不會(huì)造成燃料元件不銹鋼包殼的燒毀，但包殼溫度的快速變化與將炙熱的金屬快速浸入冷水的情況相似，會(huì)使不銹鋼包殼表面產(chǎn)生類似燒灼淬火的痕跡，在對(duì)西安脈沖堆燃料元件進(jìn)行表面觀測(cè)時(shí)，也確實(shí)發(fā)現(xiàn)了包殼表面呈藍(lán)黑色的現(xiàn)象。

圖8 D2元件燃料邊緣和包殼的溫度變化Fig.8 Temperature variation of fuel and clad for D2element

圖9示出了D2元件周圍冷卻劑通道軸向主流溫度變化情況，在開始時(shí)刻，冷卻劑軸向最高溫度并不在通道的出口處，而是位于出口下方某一位置，隨時(shí)間的增加，出口處冷卻劑溫度逐漸升高。總體分析，冷卻劑主流溫度在脈沖后上升并不顯著，這主要是由于脈沖發(fā)射持續(xù)的時(shí)間很短，實(shí)際釋放到冷卻劑中的積分能量并不高。

圖9 D2燃料元件周圍冷卻劑主流溫度隨時(shí)間的變化Fig.9 Main flow temperature variation of coolant around D2element

4.2 堆芯分析

選取D2、E2、F2、G2、H2、I2燃料元件分析反應(yīng)堆堆芯。圖10示出了6根燃料元件最高溫度的對(duì)比。D2、E2、F2元件的最大瞬時(shí)溫度要顯著高于其他元件，這是因?yàn)檫@3根元件布置于脈沖棒周圍，發(fā)射脈沖時(shí)，其功率峰因子最大。圖11示出了冷卻劑溫度沿堆芯徑向的分布，此處給出的是臨近D2、E2、F2、G2、H2、I2燃料元件的冷卻劑溫度分布圖。在堆芯內(nèi)側(cè)，由于有中央水腔，冷卻劑溫度較低；堆芯外側(cè)，因?yàn)椴贾昧舜罅康娜剂显?，脈沖發(fā)射后，通道內(nèi)的水溫變化則較大。

圖10 燃料元件最高溫度分布Fig.10 Maximum temperature distribution of fuel element

圖11 冷卻劑溫度沿堆芯徑向的分布Fig.11 Coolant temperature distribution along radius of core

5 結(jié)論

通過本文建立的脈沖堆子通道分析方法，對(duì)西安脈沖堆的脈沖運(yùn)行瞬態(tài)進(jìn)行了理論分析，結(jié)果表明：在3.5＄脈沖工況下，堆芯熱棒（D2）溫度和包殼最高溫度較大，但并未突破燃料堆芯溫度的限值，在脈沖開始階段燃料最高溫度會(huì)形成1個(gè)小的溫度峰值，但其持續(xù)時(shí)間很短，由于整個(gè)脈沖過程釋放的能量有限，不會(huì)對(duì)燃料元件造成不可接受的安全影響。

［1］ CHU P T，HOCHREITER L E，CHELEMER H，et al.THINC-Ⅳ：An improved program for thermal hydraulic analysis of rod bundle cores［M］.US：Westinghouse Electric Corporation，1973.

［2］ STEWART C W，WHEELER C L，CENA R J，et al.COBRA Ⅳ：The model and the method，BNWL 2214［R］.US：Pacific Northwest Laboratory，1997.

［3］ MERROUNA O，ALMERSB A，BARDOUNIA T，et al.Analytical benchmarks for verification of thermal-hydraulic codes based on sub-channel approach［J］.Nuclear Engineering and Design，2009，239：735-748.

［4］ CHEN Lixin，TANG Xiaobin，JIANG Xinbiao，et al.Theoretical study on boiling heat transfer in the Xi’an Pulsed Reactor［J］.Science China：Technological Sciences，2013，56（1）：137-142.

［5］ IAPWS.Release on the IAPWS international formulation 1997for the thermodynamic properties of water and steam［R］.Erlangen，Germany：IAPWS，1997.

［6］ WADE R M.Thermal hydraulic analysis of the Oregon State TRIGA Reactor using RELAP5-3D［D］.US：Oregon State University，2008.

［7］ 趙柱民，陳立新，繆正強(qiáng)，等.脈沖堆脈沖工況下堆芯功率時(shí)空分布模擬［C］∥中國核學(xué)會(huì)2009年學(xué)術(shù)年會(huì).北京：中國核學(xué)會(huì)，2009.