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

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

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

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

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

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

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

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

1.1 MOOSE 平臺(tái)簡(jiǎn)介

Multiphysics Object Oriented Simulation Environment（MOOSE）是美國(guó)愛(ài)達(dá)荷國(guó)家實(shí)驗(yàn)室開(kāi)發(fā)的面向?qū)ο蟮亩辔锢眈詈掀脚_(tái)，提供功能齊全的的底層求解器，可自動(dòng)實(shí)現(xiàn)并行計(jì)算、網(wǎng)格I/O 及多維度計(jì)算，大大簡(jiǎn)化了開(kāi)發(fā)者的工作。MOOSE 框架如圖2 所示。

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

1.2 流體傳熱模型

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

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

m˙——氦氣質(zhì)量流量；

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

其中：

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

m˙——氦氣質(zhì)量流量；

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

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

其中：

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

由公式（3）即可以得到普朗特?cái)?shù)Pr（z）：

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

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

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

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

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

1.3 固體導(dǎo)熱模型

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

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

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

1.4 流固邊界共軛傳熱模型

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

式中：q——熱流密度；

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

結(jié)合公式（13）和公式（14），公式（12）變形為：

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

1.5 PHEASANT 程序?qū)崿F(xiàn)

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

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

2 基準(zhǔn)題驗(yàn)證

2.1 MHTGR-350 基準(zhǔn)題簡(jiǎn)介

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

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

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

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

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

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

2.2 單通道計(jì)算及驗(yàn)證

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

在堆芯中，每個(gè)燃料柱均為六棱柱狀，徑向的每個(gè)冷卻劑通道周?chē)辛鶄€(gè)燃料棒，以此為最小單元規(guī)律排列。為簡(jiǎn)化計(jì)算、便于結(jié)果對(duì)比，以冷卻劑通道為中心將此單元按照面積等效原則簡(jiǎn)化成為圓形通道模型，如圖6 所示。

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

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

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

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

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

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

由于單通道模型的網(wǎng)格量較少，兩程序計(jì)算花費(fèi)的時(shí)間差距較小，因此采用組件模型進(jìn)一步驗(yàn)證該簡(jiǎn)化模型的高效性。

2.3 組件計(jì)算及驗(yàn)證

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

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

對(duì)1/12 組件模型建模畫(huà)網(wǎng)格，網(wǎng)格模型如圖10 所示，其中PHEASANT 程序采用的網(wǎng)格不包含冷卻劑通道。

在PHEASANT 和OpenFOAM 計(jì)算中，模型的三個(gè)側(cè)面和上下底面均設(shè)置為對(duì)稱(chēng)邊界條件，即絕熱邊界。冷卻劑通道與石墨接觸面的條件設(shè)置同兩程序單通道模型的計(jì)算條件，計(jì)算結(jié)果如圖11 和圖12 所示。

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

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

表6 MHTGR-350 MW 基準(zhǔn)題組件模型溫度結(jié)果Table 6 Temperature results of MHTGR-350 MW benchmark assembly model calculation

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

表7 MHTGR-350 MW 基準(zhǔn)題組件算例計(jì)算時(shí)間對(duì)比Table 7 Calculation time of the MHTGR-350 MW benchmark assembly case

PHEASANT 及OpenFOAM 均使用20 核并行計(jì)算，分別花費(fèi)3.5 min 和43 min。

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

3 結(jié)論

本文提出了“一維流體、三維固體”簡(jiǎn)化熱工模型，用于棱柱狀高溫堆堆芯熱工計(jì)算?；诙辔锢眈詈掀脚_(tái)MOOSE 開(kāi)發(fā)了PHEASANT 程序，實(shí)現(xiàn)了簡(jiǎn)化熱工模型的計(jì)算功能。

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