高温球床壁面区域有效导热系数模型优化

步珊珊，陈 波，杨光超，陈德奇，马在勇，张卢腾

(1.重庆大学 低品位能源利用技术及系统教育部重点实验室，重庆 400044；2.中广核研究院有限公司，广东 深圳 518000)

颗粒堆积球床在核能领域中有较多的应用，比如第4代核能系统中的高温气冷球床堆[1]、熔盐球床堆[2]以及聚变增殖球床包层[3]等。球床内部传热机制比较复杂，包括固体导热-表面辐射-固体导热、固体导热-气体导热-固体导热和固体导热-接触导热-固体导热等传热路径[4]，因此通常用有效导热系数来表征球床的综合传热能力。

目前球床有效导热系数模型的研究比较广泛[5-7]，其中ZBS模型[8]考虑了接触导热以及辐射，有较为广泛的应用，如Enoeda等[9]测量了钛酸锂球床、锆酸锂球床以及硅酸锂球床的有效导热系数，发现在400～800 ℃温度范围内，ZBS模型预测结果和实验测量数据吻合较好。Bu等[10]采用数值计算和实验测量研究了简单立方结构球床有效导热系数，发现ZBS模型和计算结果以及实验测量结果吻合都很好。De Beer等[11]发现ZBS模型可很好地预测球床中心区域的有效导热系数，但壁面附近区域的预测值和实验测量结果偏差较大。作者前期通过实验结果和4种不同的有效导热系数模型的对比分析也发现，ZBS模型可更好地预测球床主体区域的有效导热系数，但在壁面附近区域的预测性能需要进一步提高[12]。因此，高温球床壁面附近的有效导热系数特性需进一步研究，同时ZBS模型在球床壁面区域的预测性能也需优化。

在前期工作[10，12-14]的基础上，本文针对无序石墨球床堆的有效导热系数开展数值研究，分析无序堆积球床主体区域、近壁面区域及壁面区域有效导热系数的分布特性，针对ZBS模型在球床壁面区域的预测性能进行优化，采用前期球床实验数据及南非HTTU实验数据进行对比，分析优化后的ZBS模型在壁面区域的预测能力。

1 数值方法

1.1 计算模型

基于PFC3D软件，采用膨胀法构建无序球床堆结构：在目标孔隙率下，直径为1/2dp(dp为目标球径)的球形颗粒在指定长方体空间内随机分布，然后再将球径扩大至目标球径，通过球体颗粒间的相互挤压和碰撞运动达到平衡。为消除垂直热流方向上的壁面效应以及减少计算量，首先构建了15dp×15dp×14dp的无序堆积球床(图1a)，然后在中部区域选取6dp×6dp×14dp的球床模型(图1b)作为本文的计算模型。构建的15dp×15dp×14dp的无序堆积球床径向孔隙率分布如图2所示。模型在x方向上的孔隙率分布与De Klerk[15]提出的关联式吻合较好，平均相对误差小于5%，验证了构建的无序堆积球床结构的可靠性。

a——无序堆积球床，15dp×15dp×14dp；b——球床计算模型，6dp×6dp×14dp图1 球床堆计算模型Fig.1 Computational model of pebble bed reactor

图2 球床堆孔隙率分布及验证Fig.2 Porosity distribution and verification of pebble bed reactor

1.2 数值方法

由于球床孔隙中的氦气是静止的，计算中只考虑导热和辐射，因此只需要求解固体区域和流体区域的能量方程：

(1)

其中：ks为固体石墨材料的导热系数，取值参考文献[16]，W/(m·K)；kf为氦气的导热系数，参考NIST标准，W/(m·K)；T为温度，K；x、y、z为坐标，m。

本文计算基于ANSYS FLUENT 20.0软件平台，计算的边界条件如图3所示：左端壁面和右端壁面均为定温边界条件；另外4个壁面绝热，因此热流沿着x方向由高温壁面(左端壁面)向低温壁面(右端壁面)传递。颗粒表面之间的辐射换热采用S2S(surface-to-surface)模型计算，颗粒表面从网格单元i离开的热流qout,i计算如下：

(2)

图3 边界条件及网格划分示意图Fig.3 Boundary condition and grid system

其中：ε为发射率，本计算中壁面和球体颗粒表面均假设为漫灰体，ε均为0.8[17]；σ为Stefan-Boltzmann常数；Ti为温度，K；ρi为反射率，0.2；Fij为网格单元i和j之间的角系数，角系数的计算是基于固体颗粒表面的网格单元；qout,j为从网格单元j离开的热流。离散格式采用二阶迎风格式，当温度的残差小于10-10判定为计算收敛。当获得球床温度分布后，通过傅里叶导热定律逆求解可获得有效导热系数。本文计算工况列于表1。

表1 计算工况Table 1 Calculated case

1.3 网格独立性验证

采用FLUENT Meshing进行网格划分，采用多面体结合四面体网格的方式划分计算域，并在球体表面以及球间短圆柱接触区域进行适当加密，网格模型如图3所示。通过改变加密区域的最小尺寸，获得了3套质量较好的网格系统，其对应的最小网格尺寸分别为2.8%dp、2%dp及1%dp，对应网格的数量分别为816万、1 317万及2 654万。将上述3套网格在工况4下进行敏感性分析，计算结果列于表2，第2套网格与第3套的相对偏差仅为0.5%，可认为第2套网格达到了独立性要求，因此选择该套网格用于数值计算分析。

表2 网格独立性验证Table 2 Grid independence test

2 结果分析与讨论

2.1 局部有效导热系数分析

为计算局部有效导热系数，沿热流传递方向用13个平面将球床堆均分成14个相等的区域，每个区域的有效导热系数keff计算公式如下：

(3)

其中：q为热流密度，本计算中热量沿x方向传递，每个平面的热流密度相同，可以取高温壁面或低温壁面的热流密度，W/m2；ΔT为每个区域热流上下游平面平均温度的差值，K；Δx为每个区域的长度，1个球径。

图4 局部有效导热系数和局部固体填充率分布Fig.4 Distribution of local effective thermal conductivity and local solid fraction

局部有效导热系数和局部固体填充率(沿着x方向等分14个区域，每个区域的填充率)分布如图4所示，在球床中心区域(距离壁面4dp～11dp的区域)，填充率波动比较小，同时这个区域对应的局部有效导热系数是线性变化的，可认为这个区域是球床主体区域。在距壁面1dp之内，球床局部有效导热系数和局部填充率的变化比较剧烈，梯度最大，因此可认为是球床壁面区域；在距高温壁面1dp～4dp范围或低温壁面1dp～3dp范围内，定义为近壁面区域。在中心区域，局部固体填充率基本稳定在0.612，但局部有效导热系数均匀减小，尤其是在高温工况，这是因为局部有效导热系数分布特性不仅受到局部孔隙结构的影响，还会受到温度影响[18]。当温度较低时，导热是主要传热机制，随温度升高，石墨的导热系数减小，导热机制的贡献减小，而辐射传热机制的贡献增大。当温度超过690 K左右时，辐射逐渐成为主导传热机制[10]，因此对于工况1和2，导热是主导传热机制，主体区域的导热系数从高温壁面到低温壁面呈现上升趋势；对于工况3，靠近高温壁面的区域，辐射传热是主导传热机制；靠近低温壁面的区域，导热是主导传热机制；因此相对于工况1和2，工况3有效导热系数的增大幅度呈减小趋势。对于工况4～9，辐射传热成为主导传热机制，因此在主体区域，工况4～9有效导热系数从高温壁面到低温壁面呈现下降趋势。另外，与低温壁面区域的有效导热系数相比，工况3～9高温壁面区域的有效导热系数的增大更明显，这也是因为高温壁面区域的主导传热机制是辐射传热。

球床局部有效导热系数随温度的变化如图5所示，近壁面区域有效导热系数存在一定的分散性特点，但其随温度的变化趋势与主体区域一致；相同温度下，壁面区域有效导热系数相对于主体区域及近壁面区域明显降低，最大降低为23.8%，最小降低为16.8%，平均降低约为22%。在相同温度下，壁面区域有效导热系数的减小主要是由于固体填充率较小导致的。高温壁面区域和低温壁面区域的局部固体填充率分别为0.510和0.534，与主体区域平均填充率0.612相比，分别降低了14%和20%。壁面区域局部固体填充率小，石墨导热系数远大于氦气导热系数，因此局部有效导热系数也会更小；另一方面壁面区域局部固体填充率小，固体颗粒表面之间的辐射传热也会减弱，会进一步减小局部有效导热系数。因此，壁面区域的填充率降低会同时削弱导热和辐射传热机制，使壁面区域的有效导热系数平均下降22%。

图5 局部有效导热系数随温度的变化Fig.5 Local effective thermal conductivity vs. temperature

2.2 壁面区域有效导热系数模型优化

有效导热系数的计算结果和前期改进的ZBS模型[14]对比如图6所示，前期改进的ZBS模型对球床主体区域及近壁面区域的有效导热系数的预测性较好，其预测相对偏差大部分在15%以内，最大相对偏差仅在20%左右，且随着球床温度升高，ZBS模型的平均预测相对偏差降低；然而对于壁面区域，ZBS模型的预测值明显偏高。这说明ZBS模型未能预测壁面区域有效导热系数降低的特性。

图6 有效导热系数的计算结果和前期改进的ZBS模型结果对比Fig.6 Comparison of effective thermal conductivity calculated result and previous improved ZBS model result

图7 优化后的ZBS模型和数值计算结果对比Fig.7 Comparison of optimized ZBS model result and numerical result

根据数值计算结果(图5)，球床壁面区域有效导热系数相比近壁面区域降低22%，因此，引入修正系数Cw对改进的ZBS模型在球床壁面区域的预测值进行优化，即：

keff,wall=Cwkeff,ZBS

(4)

其中，对于球床主体区域及近壁面区域，修正系数Cw=1；对于壁面区域，修正系数Cw=0.78。优化后ZBS模型和数值结果拟合曲线的对比如图7所示，当温度小于1 600 K时，ZBS模型和计算结果在主体区域的最大相对偏差为8.7%，在壁面区域的最大相对偏差约为8.0%。需要注意的是，当温度大于1 600 K时，ZBS模型与数值计算结果的趋势有所不同。因此，本文提出的修正系数适用范围在温度小于1 600 K，与数值计算结果的相对偏差在8%以内。

2.3 优化后的壁面区域模型验证

为验证优化后的ZBS模型在壁面区域有效导热系数的预测能力，采用前期无序球床实验结果[12]和南非HTTU实验数据[11]分别进行对比，结果如图8所示。由图8a可见，与前期实验结果相比，优化后的ZBS模型在主体区域、近壁面区域及壁面区域的有效导热系数预测能力均较好，优化后的ZBS模型可预测95%以上的实验数据，其相对误差小于15%；和前期实验结果的最大相对误差为20.1%。由图8b可见，优化后的ZBS模型在主体区域的预测值与实验结果吻合较好，最大相对误差为13.1%。对于壁面及近壁面区域，ZBS模型的预测值虽仍有一定程度的偏低，但与球床有效导热系数的分布规律基本吻合。上述验证结果表明优化后的ZBS模型对球床壁面区域有效导热系数的预测能力较好。

a——与前期实验结果[12]对比；b——与南非HTTU实验数据[11]对比(82 kW)图8 优化后的ZBS模型和实验数据对比Fig.8 Comparison of optimised ZBS model and experimental data

3 结论

本文针对无序石墨球床堆有效导热系数开展了数值研究，分析了无序堆积球床主体区域、近壁面区域及壁面区域有效导热系数的分布特性，所得主要结论如下。

1) 在距壁面1dp之内，球床局部有效导热系数的梯度最大，被定义为球床壁面区域；在距离高温壁面1dp～4dp范围或低温壁面1dp～3dp范围内，定义为近壁面区域。而球床主体区域则在距离壁面大于4dp～11dp的区域。

2) 近壁面区域有效导热系数存在一定的分散性特点，但其随温度的变化趋势与主体区域一致；由于受到局部固体填充率分布的影响，壁面区域有效导热系数相对于主体区域及近壁面区域明显降低约为22%。

3) 针对前期改进后的ZBS模型，通过引入修正系数Cw来优化其在球床壁面区域的预测能力。对于球床主体区域及近壁面区域，修正系数Cw=1，对于壁面区域，修正系数Cw=0.78。通过与前期无序球床实验数据及南非HTTU实验数据的对比，验证了优化后的ZBS模型能较好地预测球床壁面区域有效导热系数。