棱柱狀高溫堆堆芯高效熱工計算模型開發(fā)及驗證

郭娟娟，童 樸，孔德泰，何 昉，黃善仿，王 侃

（1.中核戰(zhàn)略規(guī)劃研究總院有限公司，北京 100048；2.清華大學，北京 100084）

高溫氣冷堆作為GIF 公布的六種第四代反應堆候選堆型之一[1]，具有固有安全性高、發(fā)電效率高、可用作核能制氫等特點。國際上廣泛研究的高溫堆有球床式和棱柱式兩種形狀的堆芯，均使用TRISO 燃料球，采用氦氣作為冷卻劑，石墨作為慢化劑[2]。其中，由清華大學研發(fā)的球床模塊式高溫氣冷堆研究成熟，已在山東榮成建造并成功進行一次冷試，其堆芯采用球床式，燃料球由TRISO 顆粒組成，冷卻劑氦氣與燃料直接接觸；棱柱狀高溫堆的堆芯則采用石墨棱柱作為基體，如圖1 所示，基于TRISO 顆粒的圓柱狀燃料芯塊置于石墨基體的燃料通道中，氦氣不和燃料直接接觸，通過冷卻劑通道將熱量帶走。

高溫堆堆芯的熱工計算方法根據(jù)堆芯形狀有所不同，研究較為成熟的球床堆堆芯的冷卻劑和燃料球直接接觸，熱工計算通常采用多孔介質法[3]，建模時不區(qū)分燃料和冷卻劑；但是棱柱狀高溫堆堆芯的冷卻劑和燃料棒不直接接觸，橫向溫差大，多孔介質法不適用，需要對固體及氣體分別建立方程，一般采用CFD 方法計算。但是CFD 方法計算量大，計算效率低，全堆尺度的模型計算時間過長，并且流體的三維NS 方程計算時間占比較大。由于棱柱狀高溫堆堆芯的冷卻劑通道狹長，湍流強度大，單個冷卻劑通道的徑向溫差小，因此可以將冷卻劑通道的徑向溫度作平均，等效成為“一維流體”模型，而燃料和石墨基體的固體求解仍然采用三維導熱模型計算。

本文針對棱柱狀高溫堆堆芯熱工計算現(xiàn)有的問題，提出“一維流體、三維固體”熱工模型，基于國際上通用的MOOSE 平臺[4]，開發(fā)計算程序，計算MHTGR-350 基準題的單通道及組件級別算例，并與CFD 軟件OpenFOAM[5]的計算結果進行對比，驗證該模型的準確性和高效性。

1 基于MOOSE 的熱工計算模型開發(fā)

本小節(jié)詳細描述基于MOOSE 平臺開發(fā)“一維流體、三維固體”簡化熱工模型的過程，包含流體傳熱、固體導熱以及流固共軛傳熱模型。“一維流體”是指將圓形冷卻劑通道的徑向溫度作平均，得到軸向溫度隨通道高度變化的分布曲線；“三維固體”是指燃料及石墨基體采用經典的三維導熱模型求解，得到三維的固體溫度分布；而流體和固體之間的共軛傳熱用邊界條件簡化處理，利用經驗關系式求出流固接觸面的湍流換熱系數(shù)，并與流體的軸向溫度分布作為三維固體導熱求解的邊界條件。

本文的熱工模型用于棱柱狀高溫堆耦合計算，因此將程序命名為PHEASANT（Prismatic High Temperature Reactor Coupling System of Neutronics and Thermal-Fluid），包含中子物理等計算功能。

1.1 MOOSE 平臺簡介

Multiphysics Object Oriented Simulation Environment（MOOSE）是美國愛達荷國家實驗室開發(fā)的面向對象的多物理耦合平臺，提供功能齊全的的底層求解器，可自動實現(xiàn)并行計算、網格I/O 及多維度計算，大大簡化了開發(fā)者的工作。MOOSE 框架如圖2 所示。

求解器及l(fā)ibmesh 庫位于最底層，提供底層算法、網格I/O、有限元離散庫、并行計算等功能；MOOSE 框架位于中間層，提供方程的殘差余量、對角元計算等功能；位于最頂層的計算模塊或程序則需要用戶開發(fā)，用戶需要將求解問題對應的方程整理成為kernels、邊界條件、初始條件等元素，完成模型開發(fā)，實現(xiàn)最終的計算功能。

1.2 流體傳熱模型

流體的軸向溫度分布采用能量守恒方程（1）求解，

式中：cp——氦氣比熱容；

m˙——氦氣質量流量；

Q˙(z)——z高度的功率密度，這三個變量均為已知量，待求的ΔT(z)為z高度處的流體溫差，根據(jù)已知的進口流體溫度，即可得到冷卻劑通道軸向任意高度處的流體溫度分布T(z) 。湍流換熱系數(shù)采用經典的Dittus-Boelter 關系式求解，關系式的推導涉及到的參數(shù)有流體密度ρ(z) 、動力粘度μ(z) 和熱導率λ(z) ，本文將以上三個參數(shù)擬合成為T(z) 的多項式。以下為經驗關系式的推導過程：

其中：

v(z)——氦氣在高度z處的流速；

m˙——氦氣質量流量；

s——冷卻劑通道的徑向面積。

由公式（2）即可得到雷諾數(shù)Re(z)：

其中：

d——冷卻劑通道的直徑。

由公式（3）即可以得到普朗特數(shù)Pr（z）：

將公式（3）和（4）代入Dittus-Boelter 湍流換熱關系式（5）中，得到Nu(z)：

根據(jù)公式（6）可以求得換熱系數(shù)h（z）。接下來對換熱系數(shù)進行修正，利用已知的熱流密度分布q（z），求出壁面固體溫度隨高度的變化Tw（z），如公式（7）所示：

利用上式求出的壁面溫度可以得到修正系數(shù)Ct(z)：

以上的h'(z) 和T(z) 是流固共軛傳熱邊界所需要的兩個變量，作為固體導熱求解的邊界條件。

流體計算模型屬于利用已知條件和顯式公式推導變量的過程，沒有實際的數(shù)值計算，因此在MOOSE 計算模型中，公式（1）～（10）屬于輔助變量系統(tǒng)，方程中的項是Auxkernel，T(z)，h'(z)變量是Auxvariable，在數(shù)值計算開始時程序已計算出輔助變量，作為數(shù)值求解的邊界條件。

1.3 固體導熱模型

“三維固體”模型采用經典的導熱方程，如公式（11）所示。求解的固體區(qū)域包含燃料及石墨基體，本文求解單通道及三維組件級別的小尺度模型，暫不考慮輻射對傳熱的影響。

Q˙——燃料區(qū)域的功率。

在MOOSE 中，將方程（11）按照Galerkin有限元離散方法處理為弱解形式的殘差方程[6]，得到公式（12）：

1.4 流固邊界共軛傳熱模型

熱工計算中，流體和固體之間的共軛傳熱十分重要，結合1.2 節(jié)和1.3 節(jié)的推導過程，固體壁面處的熱流密度和換熱系數(shù)及流體溫度之間的關系可以表示為：

式中：q——熱流密度；

h——換熱系數(shù)h'(z) 。

結合公式（13）和公式（14），公式（12）變形為：

公式（15）即為簡化模型的核心方程，實現(xiàn)了流體和固體之間的共軛傳熱。

1.5 PHEASANT 程序實現(xiàn)

PHEASANT 程序的框架如圖3 所示，最底層的MOOSE 框架支撐程序實現(xiàn)相應計算功能，經過有限元離散的方程中的kernels、BCs 等分別歸類到PHEASANT 程序相應的塊中，增加或改變問題方程只需要在相應塊中增加或更改類即可。因此，不同方程若含有相同的kernels 或者邊界條件類型等，可以通用。

本文中，公式（15）中的kernels 和邊界條件BCs 分別在PHEASANT 相應的塊中開發(fā)好，材料的密度、導熱系數(shù)等參數(shù)在material 塊中給定，選擇Newton 求解器并行計算。

2 基準題驗證

2.1 MHTGR-350 基準題簡介

MHTGR-350 MW 基準題[7]的原型堆是MHTGR-350 MW 概念堆，由美國、韓國、德國等國家的高校或國家實驗室共同參與開發(fā)，致力于為棱柱狀高溫堆堆芯模擬計算建立相同的計算條件，從而比較不同程序之間的計算結果。該基準題熱功率350 MW，屬于棱柱狀高溫氣冷堆，如圖4 所示，氦氣作為冷卻劑，石墨作為慢化劑，每個組件都是六棱柱狀，由TRISO顆粒制成的棒狀燃料插入燃料通道中，氦氣通過冷卻劑通道將熱量帶走。

本文計算基準題的熱工穩(wěn)態(tài)算例，采用給定的功率分布及流體固體物性參數(shù)。已知的功率分布如圖5 所示，是基準題的物理臨界算例的功率分布結果[8]，該算例為平衡循環(huán)壽期末的臨界問題，仍有部分控制棒插入，因此功率曲線不完全滿足余弦分布。

計算所需的邊界條件列于表1，是基準題熱工穩(wěn)態(tài)算例設定的參數(shù)。氦氣、燃料和石墨的物性參數(shù)分別列于表2 和表3 壓力取6.39 MPa 后擬合成溫度T（單位為K）的多項式。表3 中的石墨物性參數(shù)采用全堆平均溫度所對應的值[7]。

表1 MHTGR-350 MW 基準題熱工計算主要參數(shù)[7]Table 1 Thermal-hydraulics properties in the steady-state calculation of MHTGR-350 MW benchmark

表2 MHTGR-350 MW 基準題的氦氣物性參數(shù)[7]Table 2 Helium properties of MHTGR-350 MW benchmark[7]

表3 MHTGR-350 MW 基準題燃料及石墨物性參數(shù)[7]Table 3 Fuel and graphite properties of MHTGR-350 MW benchmark[7]

2.2 單通道計算及驗證

本小節(jié)以CFD 軟件OpenFOAM 為基準，采用單通道級別算例驗證本文提出的“一維流體、三維固體”熱工模型的準確性。

在堆芯中，每個燃料柱均為六棱柱狀，徑向的每個冷卻劑通道周圍有六個燃料棒，以此為最小單元規(guī)律排列。為簡化計算、便于結果對比，以冷卻劑通道為中心將此單元按照面積等效原則簡化成為圓形通道模型，如圖6 所示。

等效之后的單通道模型以冷卻劑通道為中心，直接與中間石墨環(huán)相接觸，最外環(huán)是燃料，等效之后的軸向高度不變，仍是活性區(qū)高度793 cm，徑向尺寸列于表4 中。

表4 MHTGR-350 MW 基準題等效單通道模型尺寸Table 4 Single channel model size of MHTGR-350 MW benchmark

分別采用CFD 程序OpenFOAM 和本文的PHEASANT 程序計算此單通道算例。其中OpenFOAM 不采用簡化模型，計算完整的三維流體傳熱、三維固體導熱及流固共軛傳熱，直到滿足收斂準則。本文的PHEASANT 程序則利用開源優(yōu)勢采用本文開發(fā)的“一維流體、三維固體”簡化熱工模型。

兩個程序計算的單通道模型均采用ANSYS軟件劃分網格，由于本文的簡化熱工模型不對流體區(qū)域做數(shù)值計算，因此PHEASANT 計算的網格模型沒有流體區(qū)域，如圖 7 所示，而OpenFOAM 的計算模型有流體區(qū)域。

兩個程序計算時，除流固接觸面的邊界條件不同外，其余邊界條件、物性參數(shù)等條件完全相同。冷卻劑和燃料溫度的計算結果如圖8 所示。

圖8 中的冷卻劑及燃料軸向溫度均經過徑向平均，兩個程序的計算曲線基本吻合，冷卻劑溫度曲線的平均偏差為1.13%，最大偏差為2.30%；燃料溫度的平均偏差為0.44%，最大偏差為1.01%，擬合曲線的最大偏差均小于工程應用的5%不確定度。因此，單通道的計算結果驗證了本文提出的“一維流體、三維固體”熱工模型的準確性。

由于單通道模型的網格量較少，兩程序計算花費的時間差距較小，因此采用組件模型進一步驗證該簡化模型的高效性。

2.3 組件計算及驗證

圖9 中，紅色圓點代表燃料棒通道，藍色圓點代表氦氣冷卻劑通道，六邊形的六個綠色頂點代表可燃毒物通道，詳細的組件尺寸等參數(shù)列于表5 中。

表5 MHTGR-350 MW 基準題組件模型尺寸[7]Table 5 Size information of the MHTGR-350 MW benchmark assembly model[7]

對1/12 組件模型建模畫網格，網格模型如圖10 所示，其中PHEASANT 程序采用的網格不包含冷卻劑通道。

在PHEASANT 和OpenFOAM 計算中，模型的三個側面和上下底面均設置為對稱邊界條件，即絕熱邊界。冷卻劑通道與石墨接觸面的條件設置同兩程序單通道模型的計算條件，計算結果如圖11 和圖12 所示。

圖11 定性展示了兩程序計算的堆芯出口處的徑向溫度分布情況，除PHEASANT 沒有冷卻劑通道外，其余位置的溫度分布基本相同。圖12 中的軸向溫度經過徑向平均，兩曲線基本吻合，平均偏差為0.96%，最大偏差為1.50%。

反應堆設計重點關注平均溫度及最值溫度參數(shù)，兩程序計算得出的冷卻劑出口平均溫度、三維組件模型中燃料及石墨平均溫度、燃料最高溫度如表6 所示。

表6 MHTGR-350 MW 基準題組件模型溫度結果Table 6 Temperature results of MHTGR-350 MW benchmark assembly model calculation

根據(jù)能量守恒，兩程序計算得出的出口處冷卻劑平均溫度相等，均為687.1 ℃；兩程序計算出的燃料平均溫度、最高溫度及石墨平均溫度偏差分別為0.72%、0.28%和0.57%，最大不超過1%。對兩程序計算組件級別算例花費的時間進行比較，結果如表7 所示。

表7 MHTGR-350 MW 基準題組件算例計算時間對比Table 7 Calculation time of the MHTGR-350 MW benchmark assembly case

PHEASANT 及OpenFOAM 均使用20 核并行計算，分別花費3.5 min 和43 min。

由于PHEASANT 計算的組件模型不包含流體區(qū)域，因此網格量比OpenFOAM 少將近1/3，計算中也不需要對流體區(qū)域做數(shù)值計算；此外，PHEASANT 程序是基于MOOSE 平臺開發(fā)，計算采用Newton 算法求解方程組，相比OpenFOAM 采用的SIMPLE 算法計算效率相對更高。因此采用“一維流體、三維固體”簡化模型的PHEASANT 的計算時間比OpenFOAM節(jié)省了90%以上。

3 結論

本文提出了“一維流體、三維固體”簡化熱工模型，用于棱柱狀高溫堆堆芯熱工計算。基于多物理耦合平臺MOOSE 開發(fā)了PHEASANT 程序，實現(xiàn)了簡化熱工模型的計算功能。

采用PHEASANT 程序計算單通道及組件級別算例，并與CFD 程序OpenFOAM 進行比較。兩者計算結果的平均偏差1%，最大偏差不超過3%，且由于PHEASANT 程序所用網格量減少1/3，計算時間減少了90%以上。對比結果充分驗證了簡化熱工模型的準確性和高效性。