沈丹紅，董 博，楊 婷，傅孝良，程 旭

(1.上海交通大學 核科學與工程學院，上海 200240；2.中國核動力研究設計院 核反應堆系統(tǒng)設計技術重點實驗室，四川 成都 610213；3.國家電投集團科學技術研究院有限公司，北京 102209；4.國家能源核電軟件重點實驗室，北京 102209；5.上海電力大學 自動化工程學院，上海 200090)

在反應堆設計和安全分析中，堆芯中冷卻劑的流動和傳熱情況至關重要。燃料組件中最常用的結構是棒束結構。冷卻劑流經(jīng)棒束間的通道帶走由棒束產(chǎn)生的熱量，從而保證燃料棒四周的溫度不超過其限值。因此，增強或者改進燃料棒與冷卻劑之間的換熱能力成為反應堆熱工水力設計的重要任務。

子通道程序是目前應用較為廣泛的安全分析和熱工水力計算程序[1]。與單通道模型不同，子通道模型考慮了相鄰通道間質(zhì)量、能量和動量的橫向交換。這種考慮使熱通道的冷卻劑溫度比未考慮橫向交混時有所降低，相應的燃料元件的表面溫度與中心溫度也略有降低。因此，子通道方法提高了熱工設計的準確度和反應堆設計的經(jīng)濟性。

目前，國內(nèi)外已有大量可用于反應堆熱工水力計算的子通道程序[2-5]。本文選用COBRA-Ⅳ程序[6]對棒束中燃料棒表面溫度進行了模擬，并根據(jù)比對結果對換熱模型和空泡份額模型進行了分析。

1 實驗工況及子通道建模介紹

本文所選實驗數(shù)據(jù)來源于棒束內(nèi)流動與傳熱行為試驗[7]。試驗棒束為5×5正方形排列棒束，如圖1所示。25根加熱棒外徑D=9.5 mm，棒中心距P=12.6 mm，矩形通道截面尺寸66.1×66.1 mm。試驗段發(fā)熱長度2785 mm。自入口每隔相等距離安放一個定位格架，共6個，格架形式為普通格架，阻塞比為0.24。

圖1 棒束通道橫截面尺寸示意圖Fig.1 Cross sectional of the rod bundle

由圖1棒束通道的橫截面可以看出，棒束的排列具有對稱性，故僅取橫截面的四分之一進行子通道計算。子通道的劃分如圖2所示，共九個通道，九根燃料棒。其中，3號、6號、7號、8號為完整燃料棒的二分之一，9號為完整燃料棒的四分之一。本文中均以5號棒為研究對象進行計算值與實驗值的對比。子通道輸入卡的相關參數(shù)如表1所示。

圖2 子通道的劃分Fig.2 Cross sectional of sub-channel division

表1 子通道輸入卡的相關參數(shù)Table 1 Parameters of sub-channel input

根據(jù)實驗中的流動工況，本文對不同壓力、不同熱流密度和質(zhì)量流量的棒束流動進行了計算。表2中分別列出了單相和兩相工況下的子通道入口邊界條件的大致范圍。其中，單相工況12個，編號為D01至D12；兩相工況18個，編號為L01至L18。

表2 不同工況的入口邊界條件Table 2 Inlet boundary conditions for different working conditions

2 單相工況計算結果分析

為了較為準確地對燃料棒壁面溫度進行模擬，需要選擇合適的對流換熱模型來正確模擬燃料棒壁面與冷卻劑流體之間的換熱。由于現(xiàn)階段對于流體縱向流過棒束時換熱特性研究較少，在一般計算中，仍按照普通圓管的換熱公式進行計算，并將棒束通道的當量直徑取為特征尺寸。其中，最為常用的是Dittus-Boelter關系式(簡稱D-B公式)，如公式(1)所示，為子通道程序中所用D-B公式：

Nu=0.023·Re0.8·Pr0.4

(1)

對于棒束通道而言，不同研究者提出了關于D-B公式的修正形式。如公式所示，Weisman[8]提出了流體是水時的計算公式：

Nu=C·Re0.8·Pr1/3

(2)

式中：P——燃料棒中心距,mm；

D——燃料棒直徑,mm。

在本文的算例中，P=12.6 mm，D=9.5 mm，故C=0.031 7。

分別應用D-B公式和Weisman公式，對單相的十二個工況進行計算。每個工況中，分別計算兩個公式在子通道中的計算結果與實驗值的誤差，并求平均，做出不同工況下的計算誤差分析圖，如圖3所示。

圖3 單相各工況平均誤差分析Fig.3 Average error of single-phase operating conditions

圖3所示的十二個工況中，Weisman公式的計算誤差小于D-B公式的情況共有八個，可初步判斷選用Weisman公式得到的計算結果更好。另選取D02和D05這兩個工況進行分析，如圖4所示。

圖4 單項工況計算結果對比圖Fig.4 Comparison of calculation results in single phase conditions

由圖3的誤差分析可知，D02和D05這兩個工況中，兩個公式所得計算誤差均較其他工況偏大。D02工況中，兩個工況均有部分點與實驗工況較為吻合，而總體上看D-B公式的趨勢與實驗值更接近，故該工況下D-B公式平均計算誤差更小。在D05工況中，兩個公式計算所得結果與實驗值相比均較高，而使用Weisman公式所得結果更接近實驗值。與其誤差結果相似的D09工況中，也是如此，在計算結果偏離實驗值較遠時，使用Weisman公式可以得到更為接近實驗值的結果。

綜合以上分析，在本文進行的計算范圍內(nèi)，與D-B公式相比，選用Weisman公式進行計算能夠在整體上獲得更好的計算結果。

3 兩相工況計算結果分析

3.1 兩相換熱模型

在兩相流動中，冷卻劑與壁面之間的換熱情況非常復雜。在進行模擬計算時，不僅要判斷冷卻劑的換熱工況，還要給出各個工況下合適的換熱模型。在較早版本的COBRA程序中，兩相換熱系數(shù)完全由用戶提供[9,10]。但在目前的COBRA程序中，有一個來源于RELAP4的公式包可供選擇[11]。

選用該公式包進行計算時，程序會根據(jù)壁溫、空泡份額等條件判斷冷卻劑的換熱工況，從而選擇合適的換熱公式。表3列出了該公式包所劃分的七個流動區(qū)域及相應的換熱公式名稱。

表3 COBRA中的兩相換熱公式[11]Table 3 Two-phase heat transfer equations in COBRA

以工況L01為例，如圖5所示，實線為選用RELAP4公式包所算得的壁面溫度變化。結合該公式包的判斷標準，對該曲線解讀如下所示。

圖5 RELAP4公式包計算所得壁溫變化對比圖Fig.5 Comparison of wall temperature calculated by RELAP4 formula package

(1)在壁面溫度達到飽和溫度之前，程序判定為強制對流階段，默認采用D-B公式進行計算。通過對單相結果的分析，本文中均采用Weisman公式進行計算，故此階段所用公式為單相時的換熱公式，即本文所選Weisman公式。

(2)壁面溫度達到飽和溫度之后，換熱公式采用Thom提出的換熱模型[12]，如公式(3)所示：

(3)

式中：ΔTsat——壁面溫度與流體飽和溫度的差；

q——熱流密度;

p——系統(tǒng)壓力。

公式(3)計算所得溫度差為華氏度單位。由該公式可以看出，此時壁面與流體間的換熱僅與熱流密度和壓力有關，在此工況下為定值。由計算結果也可看出，此時壁面溫度不變化，該公式無法模擬出定位格架的效應。

(3)壁面溫度變化的后半段，有所升高，此時不再使用Thom換熱模型。當空泡份額大于0.8時，換熱模型改換為Schrock and Grossman模型[13,14]，如公式(4)至公式(6)所示。

(4)

(5)

(6)

該模型雖可以模擬定位格架效應，但計算結果與實驗數(shù)據(jù)不符，說明對于流動區(qū)域的判斷有誤。而該判斷準則是依據(jù)空泡份額的大小來進行流動區(qū)域的判定的。因此，對于空泡份額的相關計算模型需要進行改進。

3.2 空泡份額模型

空泡份額α是兩相流的基本參數(shù)之一，由于它的測量比較困難，因此必須借助空泡份額的關系式進行計算。在上一節(jié)所示的計算中，選用的空泡份額模型為均相模型，其結果顯示該模型并不適用于此工況的計算。為了找到合適的空泡份額模型，用COBRA程序中的幾個模型分別進行計算，并將模型的具體形式列如表4所示。

表4 空泡份額模型公式[15]Table 4 Correlations of void fraction model

表4所列公式中，vG、vL分別為汽相和液相的流速，x為質(zhì)量含氣率(0≤x≤1)，S為滑速比(S>1)。通過對公式的簡單分析可知，同一工況下，用均相模型計算所得空泡份額數(shù)值最大，其余兩個模型計算結果較之稍小。

仍以工況L01為例，用不同模型計算所得空泡份額變化如圖6(a)所示，所得壁溫變化如圖6(b)所示。

圖6 不同空泡份額模型計算結果對比Fig.6 Comparison of different void fraction models

由圖6(a)可看出，與從公式進行分析得到的結論相同，均相模型得到的計算結果數(shù)值最大，Slip模型所得結果隨滑速比S的增大而減小。S=3.0時的Slip模型和Modified Armand模型所得結果較小，然而兩者所得空泡份額數(shù)值變化趨勢不同，A05所得結果最終未超過0.8，而其余幾組結果均在不同位置空泡份額的數(shù)值超過了0.8。根據(jù)RELAP4公式包的判斷準則，當空泡份額超過0.8時，不再使用Thom換熱模型，而選擇Schrock and Grossman模型。故在圖中所示的計算所得壁溫有著不同的變化趨勢。其中，與實驗值吻合最好的計算結果為A04，A05，相對應的空泡份額模型分別為Slip模型和Modified Armand模型。由于使用Slip模型時，在不同工況下選用精確的滑速比是較為困難的，故本文最終選用Modified Armand模型進行計算。

3.3 兩相計算結果

對兩相的十八個工況進行計算，并選擇其中三個工況為一組進行分析。如圖7所示，分別為工況L15、L01和L08。這三個工況質(zhì)量流率和熱流密度一樣，壓力分別為1.5 MPa、2.0 MPa和4.0 MPa。

圖7 兩相工況計算結果對比圖Fig.7 Comparison of calculation results in two phase conditions

由圖7可以看出，在本文所涉及的兩相工況下，壁溫上升到一定溫度后保持不變。計算所得壁溫較實驗值更早的出現(xiàn)這種變化。這是由于在單相部分，計算值均略高于實驗值。這使得計算所得壁溫較早的達到飽和溫度，換熱模型由單相換熱模型變?yōu)門hom換熱模型，同一工況下，由該模型計算所得壁溫為定值。

對計算結果進行誤差分析，上述L15、L01、L08工況下的平均誤差分別為0.79%、0.89%、0.86%，可見Thom換熱模型對于壓力這一參數(shù)的變化并不敏感。由于兩相工況下的誤差一部分來源于單相換熱模型，一部分來源于Thom換熱模型，且單相換熱模型的誤差更大，因此平均誤差的分析對于Thom換熱模型而言意義不大。Thom模型主要影響了兩相階段計算值與實驗值的偏差，所有工況所得計算結果與實驗值的溫度差均不超過10 ℃，誤差小于5%，計算結果及對比圖不在此一一列舉。

綜上，兩相工況下選用RELAP4公式包進行換熱模型的選擇?？张莘蓊~作為公式包進行狀態(tài)判斷的重要參數(shù)之一，需要選擇合適的模型對其進行計算。經(jīng)過對比，本文選擇Modified Armand模型計算空泡份額。兩相換熱模型中的Thom換熱模型能夠得到較好的換熱結果。但單相換熱模型對兩相計算仍有較大影響，需要進一步改進單相換熱模型。

4 結論

本文利用子通道程序COBRA對燃料棒壁面溫度進行模擬，通過對計算結果與實驗數(shù)據(jù)進行對比，對不同模型的適用性有了一定的認識，結論如下：

(1)單相對流換熱公式D-B公式和Weisman公式相比較下，Weisman公式在整體上能夠獲得較好的計算結果;

(2)兩相工況選用RELAP4公式包進行計算，計算結果的準確性取決于公式包對流動工況判斷的準確性及各個換熱模型的準確性。選用Modified Armand模型計算空泡份額，使得對流動工況的判斷較為準確;

(3)在兩相計算時，單相換熱模型對計算結果存在較大影響，單相換熱模型仍需要進一步改進。