球床式高温气冷堆是第四代核反应堆类型之一。在球床反应堆中，氦气作为冷却剂，核燃料嵌在石墨慢化剂中并可以维持很高的温度。这种结构可以保证燃料不会熔化，使反应堆具有良好的固有安全性^［1-4］。该堆型具有结构简单、燃料元件适合批量生产以及燃料球装卸方便等特点。

目前，很多学者对球床中的流动和换热特点进行了研究。孟现柯等^［5］用由碳钢球堆积成球床，以蒸馏水为工质，采用电磁感应加热方式研究了球床通道内部的换热特性；Jia等^［6］用黑色玻璃球堆积成球床，以水为工质研究了不同床型对内部工质流动情况的影响。考虑到实验成本和复杂性，多数实验结果很难为球床反应堆提供可靠的参考依据。Lee等^［7］于2007年采用 ${{k}} - \varepsilon $ 湍流模型比较了不同球间间隙对流场分布的影响；2014年，Ferng等^［8］用雷诺应力湍流模型、非结构网格对体心立方（BCC）人工间隙和点接触两种接触方式的共轭模型的模拟结果进行了分析比较；2012—2015年，Shams等^［9-11］采用q−DNS方法、多面体网格对面心立方（FCC）人工间隙模型进行了分析，发现在两球之间的狭窄区域存在射流现象及速度场、温度场分布的不对称性；2017年，蒋旭等^［12］采用大涡模拟（LES）方法、结构网格研究了BCC堆积模型的流场和温度场，通过与 ${{k - }}\varepsilon $ 湍流模型计算结果对比发现，采用LES方法能够更好地捕捉流场中的涡结构，瞬时温度变化更具参考价值。采用数值方法可以更容易获得更多的流场、温度场的信息，但多数学者采用数值方法时均只考虑球床在强制对流情况下的流动传热特性，而忽略了热浮力的影响。Huang等^［13］对矩形腔内湍流混合对流特性进行了实验和数值研究，分析了浮力的影响，通过与实验结果进行对比，确定了LES方法在混合对流计算中的精度；Guardo等^［14-15］采用数值方法，对随机球床内超临界二氧化碳在高压下的混合对流传质进行了分析，研究表明，对于层流状态下的超临界流体，可以直观地看到流量和流动方向对传质量的影响，在助流状态下，会得到更大的传质量。研究表明，混合对流中浮力的影响不可忽略，但对球床内流动和传热的研究还不够充分。

本文将采用LES方法、结构网格研究热浮力对球床反应堆内FCC面接触单元的燃料球表面及局部流动与传热的影响。

1 几何模型与数值方法

球的规则排列方式主要有3种^［16］：简单立方（SC），体心立方（BCC），面心立方（FCC）。球床内部燃料球堆积结构可以看作是3种方式的随机组合。两球的接触点一般有3种处理方式：人工间隙、点接触、面接触。3种规则排列方式和接触点处理方式如图1所示。由于3种排列方式中FCC结构空隙率最小，而实际球床中燃料球的无规则堆积必然使球床的空隙率最小化，因此本文选择燃料球的FCC排列方式，又考虑到两球在接触位置会存在一定程度的变形，对接触点的处理方式选择了面接触，由此得到了一个FCC单元。考虑到计算资源有限及为避免受进、出口效应影响，本文中采用4个FCC单元堆积模型，几何模型如图2所示。主要参数参考Shams等^［11］中的值。燃料球球面材料为石墨，直径为60 mm，接触变形部分是燃料球直径的1%，内部冷却工质为氦气。边界条件及参数、物性参数分别如表1、2所示。

考虑热浮力时，密度采用Boussinesq近似，计算 ${G_{\rm{r}}}/{Re}^2={g^*}\;{\alpha _{\rm{v}}}\;\Delta tl/{v^2}$ = 0.024，其中： ${g^*}$ 为本文中采用的重力加速度， ${g^*}={{g}}{\alpha _{\rm{v}}}$ （ ${T_{\max }} - $ $ {T_{\min }}$ ），T_max和T_min分别为； ${\alpha _{\rm{v}}}$ 为体积膨胀系数； $\Delta t$ 为进、出口温差； $l$ 为特征长度； $v$ 为流体进口速度；G_r为格拉晓夫数；Re为雷诺数。复合对流中，采用 ${G_{\rm{r}}}/{Re}^2$ 作为判断自然对流影响程度的判据。一般认为，当 ${G_{\rm{r}}}/{Re}^2 \;\geqslant\; 0.01$ 时，自然对流的影响不能忽略。本文中重力方向与流动方向相同，热浮力方向与流动方向相反。

为了能够捕捉更多的瞬时流场信息，采用LES方法进行研究。亚格子模型采用动态Smagorinsky−Lily模型，其优点在于具有一定的自适应性，得出的结果更可靠。

2 网格适用性验证

为了提高计算结果的准确性，对所采用的FCC面接触模型进行结构网格创建，采用ICEM软件将流场区域划分成多个六面体几何结构并建立对应的六面体块，最后通过将几何体与对应的六面体块进行关联，成功完成了结构网格的划分，并对近壁面网格进行了加密。FCC面接触单元网格如图3所示。

为了验证近壁面网格加密方法的适用性，对Shams等^［12］采用的FCC人工间隙模型，采用相同的结构网格划分方法进行全结构化网格的划分。在相同的边界条件和进、出口参数条件下，采用LES方法进行计算，获取了在垂直来流方向上近球面到中心流场区域的时均速度分布。图4为近壁面加密网格适用性验证结果。4（a）为提取数据位置。图4（b）为LES和q−DNS方法结果对比，其中，横坐标表示的位置与纵坐标表示的时均速度都进行了无量纲化，并与Shams等^［14］的q−DNS数据进行对比； $\overline U $ 为时均速度， ${\overline U _{\max }}$ 为最大时均速度；（ $Z - {Z_{{\rm{mid}}}}$ ）为在中心流场区域所提取线上的数据点Z相对中间位置Z_mid的坐标；D为燃料球直径。结果显示，与近球面处计算结果高度一致，由此说明近壁面的网格加密方法适用于该球床模型的数值计算。

为验证网格无关性，分别对网格数为720万、503万、325万的模型进行了计算，并整理对比了进、出口速度，温度差。以720万网格下的计算结果为标准，计算出另外两种网格数下得到的数据误差。网格无关性验证结果如表3所示，结果显示3种网格数的计算结果无明显差异，这说明325万网格已可满足计算要求，网格数的增加对计算结果无明显影响。由于采用了LES方法，为了能更好地捕捉内部流场信息，本文中用于计算的网格数取为503万。

3 计算结果与分析

在所研究的流场区域沿流动方向取相同的截面，获取截面上的平均速度，各层截面位置及平均速度、温度变化如图5所示。结果显示，在进、出口处速度变化较大，存在明显的进、出口效应，扰动较为剧烈。进口处受影响的程度相对较大，并持续到第3层截面之后。从第4层截面开始，流速呈现对称性波动，表明流动已基本稳定。可以发现，自第4层截面之后，热浮力影响效果明显，受热浮力影响的流体流速相对较低，说明热浮力的存在对流速起到抑制作用。第2层截面之后，每层温度均有升高，相邻两层温差ΔT为7 K。为了进一步研究流场内部的流动和传热情况，通过以上的分析，选取第5层和第7层截面之间的FCC单元作为流场研究区域，选取3号中心球分析球面换热情况。

3.1 热浮力对速度和温度场的影响

在所研究的FCC单元内选取截面1、2、3，流场区域各截面位置及截面上流线图如图6所示。获取截面1、3上的时均速度和时均温度分布，结果如图7、8所示。可以发现热浮力对局部流场速度与温度分布影响较为明显。在截面1的速度场中，受热浮力影响，低速区域增大，尤其在燃料球底部和顶部对应的尾流区和滞止区位置，低速区域明显增大。温度场中，低速区域的增大使高温区域随之增大，且在中心区域，尾流区温度较高。截面3与接触面相切，对比速度分布，在尾流区和滞止区同样可以发现与截面1中相同的现象。在接触面处热浮力存在时，流体速度会在发生分离的位置略微增加，之后的速度都会相对减小。截面上的平均温度相对较高，接触面处较为明显，高温区域延伸较长。这再次说明热浮力的存在对流速起到抑制作用。

为了解速度和温度在球面和中心流场区域之间的变化情况，在截面1上中心区域沿流动方向和垂直流动方向各取一条截线，如图9所示。提取线上的时均速度和时均温度分布，并与不考虑热浮力的分布进行了对比，结果如图10所示，其中：V_b 、T_b分别为考虑热浮力时的时均速度和时均温度；V_n 、T_n分别为不考虑热浮力时的时均速度和时均温度。分析截线1上提取到的时均速度和温度分布（相对于中心点）：沿X正方向，在近球面位置热浮力对速度和温度的影响并不明显；远离球面后，热浮力的影响开始增强。在热浮力影响下，相对中心点的对称位置处较低的速度被提高，较高的速度被降低，最终减少了速度分布的不对称性，也同样减少了温度分布的不对称性，截线上的温度有所提高，考虑到速度差的减小，降低了湍流强度，削弱了换热。相对于不考虑热浮力的情况，靠近球面位置，速度和温度的相对变化较大（速度最大降低70%，相对层间温差，温度最大增加43%）。远离球面后速度和温度相对变化趋于平缓，在中心流场区域速度相对变化范围在10%以内，温度平均增加约20%。

分析截线2上提取到的时均速度分布，沿Y负方向（流动方向），不考虑热浮力时，由上一个中心球底部到下一个中心球顶部，速度先快速增大，后缓慢增大，在中间偏下的位置达到峰值，后缓慢减少，至靠近球面位置急剧降低最后接近滞止，速度受滞止区和尾流区影响明显。在考虑热浮力后，由速度分布可以明显发现，相对中心点呈良好的对称性，截线上速度相对减小，峰值位置更靠近中心点。靠近球面处速度剧变的情况有所缓和，从一定程度上减少了由于速度变化大引起的湍流强度。从时均温度分布可以发现，考虑热浮力后，高低温差减小，中心区域截线上平均温度有所增加。靠近球面位置，时均速度和温度的相对变化较大，速度几乎滞止（相对层间温差，温度最大增加50%）。靠近中心流场区域，速度和温度的相对变化有所减小，温度最少升高约20%。

3.2 热浮力对燃料球表面对流换热的影响

图11为3号中心球表面温度分布。由图中可见，热浮力的存在使球面部分位置温度相对升高，尤其在燃料球的顶部、侧面的4个接触面中间部分以及底部较为明显。本文中定义努塞尔数Nu = $hD/\lambda $ ， $h = q/(T - {T_{\rm{in}}})$ ，其中： $h$ 为对流换热系数； $\lambda $ 为氦气热导率； $q$ 为热流密度； $T$ 为球面上不同位置温度； ${T_{\rm{in}}}$ 为氦气进口温度。球面Nu分布如图12所示。对比发现，热浮力存在时，燃料球顶部及接触面附近的Nu明显降低，高Nu区域明显减小，且差异较明显的位置分布在燃料球的上半球面。由时均速度分布情况，顶部附近流场速度相对减小，减弱了对流换热，使局部Nu降低。由图6截面1、2、3上流线图可以发现，流场中的涡多集中在接触面后沿及燃料球底部的尾流区域，扰动较为剧烈。这说明热浮力的存在削弱了换热强度，而涡的存在增强了局部扰动，削弱了热浮力的影响。

为进一步了解顶部和接触面附近的换热情况，在球面上3个特殊位置取3条截线，截线3、4、5，取线位置如图13所示。球面3条截线上Nu分布如图14所示。通过对比每条截线上的Nu分布可以发现，球面顶部换热最好，接触面附近由于接触面的阻碍作用Nu较小，换热较差。在考虑热浮力的情况下，球面相同位置的Nu会相对减少，最大差值约为球面平均Nu的6%，即热浮力的存在对对流换热起到抑制作用。沿流动方向，位于球面不同位置，热浮力影响的程度也有所不同。由图中可以看出，燃料球的上半部分受热浮力的影响较大，而下半部分由于流体扰动较强受热浮力影响较小，这与分析球面Nu分布的结果一致。

4 结　论

采用LES方法，在考虑热浮力的情况下，对球床反应堆FCC面接触结构内的流动和传热进行了研究。通过对比分析整个球床各层的平均速度变化，局部流场中速度、温度场的分布及中心燃料球表面的温度、Nu分布，得到热浮力对球床反应堆局部流动和传热影响的几点规律：

（1）低雷诺数下热浮力的影响不能忽略，热浮力的存在对球床中各层平均速度变化起到阻碍作用，会抑制速度产生较大的波动。

（2）在局部流场区域，时均速度、时均温度受热浮力影响明显。靠近球面位置，相对时均速度和相对时均温度变化较大。在中心区域，热浮力的存在会减小流场中速度分布的不对称性，速度最大降低约10%，时均温度最少升高约20%。

（3）相对来流方向，燃料球上半部分，尤其在燃料球顶部及接触面附近，受热浮力影响Nu减小，最大降低6%，削弱了换热。流体流过接触面接近球底部，涡结构增加，扰动剧烈，削弱了热浮力影响。