脈沖反應堆脈沖運行瞬態的子通道分析

陳立新，朱 磊，馬騰躍，江新標，陳 達

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

脈沖反應堆是一種以鈾氫鋯為燃料的池式研究堆，由于鈾氫鋯材料具有較大的瞬發負溫度系數，因此該類型反應堆既可穩態運行，也可脈沖運行。對于脈沖運行工況，反應堆的功率變化快，脈沖功率峰的寬度達ms量級，功率變化的范圍大，脈沖前的穩態功率與脈沖峰值功率相差6～7個量級，這種快速的功率變化對堆芯的傳熱安全提出了更高要求。準確計算脈沖運行狀態下堆芯的熱工參數是評價脈沖反應堆瞬態安全的重要內容。

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

1 子通道物理模型

1.1 子通道模型的數學方程

一般，子通道分析常采用均勻平衡模型，即假設氣、液兩相的相對速度為零，且具有相同的壓力與溫度。基于以上假設，可將冷卻劑看作氣、液兩相的混合物而進行統一處理，其連續性方程及能量和動量方程［3］如下。

連續性方程為：

能量守恒方程為：

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

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

1.2 子通道方程的數值解法

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

對于離散化方程守恒方程，在每個時間步長Δτ內按軸向分層逐層對各通道進行迭代求解，得到本層的所有參數后，進入下一層的求解，其中，能量方程需同燃料元件導熱方程聯合求解。

圖1 子通道網格劃分Fig.1 Mesh generation of sub-channel

2 輔助模型

2.1 功率計算模型

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

2.2 燃料元件導熱模型

對于燃料元件的導熱，采用粗棒非穩態一維導熱模型，該模型可用于模擬燃料元件穩態和瞬態導熱問題，其通用控制方程為：

式中：r為與熱量傳遞方向平行的坐標；F（r）為與導熱面積有關的因子；Q 為內熱源；λ為導熱系數；T 為溫度。對于式（4），可將其在時間和空間上進行離散，采用差分法進行求解。

2.3 熱傳輸模型

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

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

2.4 水和水蒸氣物性的計算模型

選用水和水蒸氣性質國際聯合會（IAPWS）發布的IAPWS-IF97公式［5］計算水和水蒸氣的物性。IAPWS-IF97公式將水和蒸氣按性質的不同分成5個區：1區為常規水區，2區為過熱蒸氣區，3區為臨界區，4區為濕蒸氣區，5區為低壓高溫區。本文主要用到1、2 區的數據。IAPWS-IF97 適 用 范 圍 為273.15 K ≤T ≤2 273.15K，壓力p≤100MPa。每個分區的物性參數計算公式參見文獻［5］。

3 計算程序的編制及驗證

3.1 程序編制

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

3.2 程序驗證

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

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

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

4 脈沖堆瞬態特性分析

為準確模擬反應堆在脈沖工況下堆芯熱工參數的時間-空間變化，采用三維時空動力學程序XAPR-HENKO 對脈沖堆在引入3.5＄脈沖時的堆芯三維功率分布進行計算［7］，脈沖工況的動力學參數列于表1。

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

表1 脈沖工況的動力學參數Table 1 Kinetics parameters in pulse operation

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

4.1 單棒計算

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

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

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

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

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

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

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

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

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

4.2 堆芯分析

選取D2、E2、F2、G2、H2、I2燃料元件分析反應堆堆芯。圖10示出了6根燃料元件最高溫度的對比。D2、E2、F2元件的最大瞬時溫度要顯著高于其他元件，這是因為這3根元件布置于脈沖棒周圍，發射脈沖時，其功率峰因子最大。圖11示出了冷卻劑溫度沿堆芯徑向的分布，此處給出的是臨近D2、E2、F2、G2、H2、I2燃料元件的冷卻劑溫度分布圖。在堆芯內側，由于有中央水腔，冷卻劑溫度較低；堆芯外側，因為布置了大量的燃料元件，脈沖發射后，通道內的水溫變化則較大。

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

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

5 結論

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

［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］ 趙柱民，陳立新，繆正強，等.脈沖堆脈沖工況下堆芯功率時空分布模擬［C］∥中國核學會2009年學術年會.北京：中國核學會，2009.