高温气冷堆(HTR)^［1-3］是第四代核反应堆，以石墨为慢化剂包覆天然铀组成球形燃料元件，氦气流通球床进行冷却，从而实现反应堆的固有安全性。固有安全性可以保证燃料不熔化，但局部高温热点依然可能破坏燃料球。因而，为了保证燃料球的完整性，寻找球床反应堆内部高温区域是十分必要的。传统实验方法在寻找堆积球内部高温点时比较困难，所以越来越多研究者开始采用数值分析方法探究球床内部温度问题。

实际球床反应堆中燃料球随机堆积，整体结构非常复杂，因而很多研究都对球床问题做了简化，并给出了四种规则排列方式：简单立方(SC)、体心立方(BCC)、面心立方(FCC)和体心棱柱(BCP)^［4］。Hassan^［5］于2007年用大涡模拟(LES)方法对BCC单层结构建立了非结构网格并进行了分析。Ferng等^［6］于2013年采用雷诺应力模型(RSM)对BCC与FCC单层结构进行了研究。Shams等^［7-8］于2013年采用近似直接数值求解(q−DNS)方法对FCC单层结构的流场温度场进行求解。蒋旭等^［9］于2017年在研究BCC结构时比较了k-ε模型与LES方法的差异。姚强等^［10］于2019年对8层排列FCC球阵建立结构化网格，并采用LES方法进行数值求解。吴浩等^［11］于2019年采用DEM模型对HTR−10球床气固传热过程进行了模拟计算，其结果基本与经验模型一致。张双宝等^［12］于2020年建立耦合模型模拟稳态流动换热过程，并将CFD结果与Thermix计算结果作比较，事实证明CFD模型精准可靠。孙世妍等^［13］于2020年研究了不同工况下HTR−10堆芯温度场的分布，以确定最高温度变化与冷却剂的关系并寻找高温极限。结合众多研究成果可以发现，规则排列球阵中，有关BCP结构的研究较少；RANS湍流模型由于捕捉瞬态流场的不足和对流动分离的低敏感性，在研究多曲面流道的球床问题上并不适用；使用结构网格在流道狭窄处可以得到更精细的结果，适合复杂的球床问题；单层结构由于整体尺度小，受出入口效应影响较大，可能最终会影响温度分布，使用多层结构并选取中间层进行研究可以较好地避免这一问题。

综上所述，本文选择四种规则排列方式中空隙率最小、排列最紧密的体心棱柱(BCP)结构，燃料球间接触位置以面接触方式处理，湍流模型选择LES方法，对8层燃料球阵建立结构网格求解温度场与速度场，寻找流场内高温区域，分析高温区域产生的原因。同时，充分发挥LES方法的优越性，寻找并分析球床内部不同位置涡对流动换热的影响。

1 计算模型及方法 1.1 计算模型

在规则球阵数值计算研究中，出入口效应的影响不可忽略，单层结构若不经过特殊设置其结果必然存在较大偏差。建立多层结构，然后选取其中间层的燃料球进行研究是一个比较好的选择。入口效应最多影响到第3层，出口效应的影响相对不明显，第4层结构受出、入口效应的影响已经不大，而采取8层结构的球阵并研究第4、5层球则能完全避免出、入口效应。本文采用的模型为面接触体心棱柱(BCP)排列，如图1所示，图中由左向右分为8层，左端为速度入口，右端为压力出口，氦气作为冷却气体从左端通入，不考虑重力。燃料球直径为60 mm，接触面设为直径的1%。燃料球面材料为石墨，且球面设为热源8 317 W·m⁻²。燃料球阵四周设为对称面，不考虑壁面效应带来的影响。物性参数参考文献［14］中的数据设置，具体如表1所示。

在选择湍流模型时，为了准确捕捉流场中的涡结构，本文采用大涡模拟(LES)方法。LES方法通过一个滤波函数对大涡直接求解，对小涡过滤，并建立亚格子应力(SGS)模型。在众多SGS模型中，本文选择了自适应更好的动态Smagorinsky−Lilly模型^［15］。网格方面，本文使用Icem软件，先寻找体心棱柱结构的最小几何单元(1/32)，然后将最小几何单元划分成11个六面体，通过点、线、面关联构建最小几何单元的结构化网格，并最终完成完整的面接触体心棱柱球阵的结构网格。BCP面接触结构网格如图2所示。最终达成单层体心棱柱球阵131万网格、8层1 050万网格以满足大涡模拟LES方法的计算要求。

1.2 计算验证

为了验证网格与计算方法的适用性，本文选取与文献[16]中同样的中心球之间截线位置，并对距离与时均速度进行了无量纲处理，对比结果如图3所示。由于模型和计算方法的差异，截线中心区域的结果差别明显，而且LES速度结果相比q−DNS速度结果更加对称。除此之外，两结果时均速度趋势一致，在高速区域和近接触面(点)低速区域结果基本一致，由此可以说明本文所采用的计算方法与网格数量可以满足计算要求。

2 结果与讨论

考虑到出、入口效应的影响，本文选择了中心位置第4层燃料球及其接触燃料球作为主要的研究对象。同时，为了便于观察各接触面附近区域的速度场温度场等细节，在体心棱柱球阵内部流域的3个位置截取截面1～3，并在截面1上取线1～5，具体位置如图4(a)、(b)所示。

图5为体心棱柱单层结构的示意图。根据接触面位置的不同可以将接触面分为：体心棱柱顶点燃料球相互接触形成的接触面1、体心棱柱同层中心燃料球相互接触形成的接触面2、体心棱柱顶部燃料球与中心燃料球相互接触形成的接触面3以及体心棱柱底部燃料球与中心燃料球相互接触形成的接触面4。

2.1 BCP结构面接触与人工间隙处理方式对比

在处理规则排列球阵接触面的问题上，人为拉开燃料球产生间隙来简化问题的处理方式比较常见。本文将对面接触与人工间隙处理方式作简单对比分析。为了便于观察，本文采用球坐标系，Φ代表天顶角，θ代表方位角。在球面取θ=0°和θ=45°两处特殊位置曲线，如图6所示。

图7为球面温度对比。人工间隙处理方式由于不存在接触，整体温度比面接触方式的低5～10 K左右。面接触结构中由于接触面的存在，曲线并不连续，其中θ=0°曲线经过接触面3，θ=45°曲线经过接触面2、4。随着氦气的流动，Φ值由正变负，在接触面曲线不连续的位置，Φ值较大的位置称为接触面前沿，Φ值较小的位置称为接触面后沿。不管是面接触还是人工间隙处理方式，原接触面位置必然会产生温度升高的现象，而面接触方式的温差更为明显。接触面前沿过渡到接触面后沿时温度骤升，且接触面后沿的温度远大于同位置人工间隙处理方式的温度。θ=0°这条曲线显示氦气流经接触面3时，对应的高温区域基本一致；而θ=45°的曲线则显示出接触面2、4处采用人工间隙处理方式时的温度趋势相比于面接触的出现了延迟。当面接触方式的温度由最高点下降后，人工间隙处理方式的温度刚刚到达最高温度，由此可见不同的接触面影响各不相同。

总体来看，采用人工间隙处理方式的温度分布更为均匀，但其高温区出现的位置与实际高温区域并不相符；面接触处理方式是两球真实接触得出的结果，温度的不均匀代表流场的复杂性。因此，面接触方式更为合理。

2.2 BCP结构流动传热特性

BCP结构内部流场由于各种接触面的存在十分复杂，本小节将从流域进行分析。

图8、9分别为流域内各处截面时均与瞬时温度和速度分布。观察温度结果可以发现，接触面附近存在温度高于其他流域温度的区域，这些高温区域集中于接触面后沿。与时均结果相比，瞬时温度分布更加混乱，呈现出一种不对称性，局部位置瞬时温度明显高于时均值。观察速度结果可以发现，在中心球顶部与底部存在低速区域，顶部区域为滞止区，底部区域为尾流区。滞止区实际上流道逐渐变宽，氦气沿着球面继续流动，顶部只流进少量氦气并附着停留；而尾流区流道逐渐收窄，同时因为接近接触面后沿，流动出现了分离，分离出来的氦气与下一中心燃料球顶部滞止区留存的氦气相互作用。这种相互作用导致氦气的再附着，因此尾流区低速面积比滞止区的更大。对比速度场与温度场可以看到，接触面位置处速度低且始终对应高温，可见流道狭窄小的接触面流动换热欠佳；而滞止区并未出现明显的高温区，可见速度分布不是影响温度分布的唯一因素。

图10为截面1取线位置时均速度。由图中可以观察到，在接触面附近存在非常明显的低速区域，而两中心燃料球之间区域接近尾流区和滞止区，速度也较小；高速区域则代表未受接触面影响的氦气脱离接触面、尾流区、滞止区以较快的速度进入下一层BCP结构。线1～5位置处的时均速度变化表明在从接触面前沿向后沿过渡时氦气速度不断下降，此时流动逐渐产生分离，这也是速度变小的根本原因。从前沿向后沿过渡时，接触面附近会产生一个较小的速度波动且随着来流方向越来越明显，这个速度波动代表接触面后沿产生了漩涡。

2.3 BCP结构中心球涡分析

为了判断涡的存在，本文将采用Q准则^［17］和涡量^［18］进行对比分析，并选取更合适的判据。Q准则由Hunt等^［17］于1988年提出，定义为

$ \qquad Q=-\dfrac{1}{2}({S}_{ij}{S}_{ij}-{\mathrm{\Omega }}_{ij}{\mathrm{\Omega }}_{ij}) $

(1)

$\qquad {S}_{ij}=\dfrac{1}{2}\left(\dfrac{\partial {u}_{i}}{\partial {x}_{j}} + \dfrac{\partial {u}_{j}}{\partial {x}_{i}}\right) $

(2)

$\qquad {\varOmega }_{ij}=\dfrac{1}{2}\left(\dfrac{\partial {u}_{i}}{\partial {x}_{j}}-\dfrac{\partial {u}_{j}}{\partial {x}_{i}}\right) $

(3)

式中： $ {u}_{i} $ 和 $ {u}_{j} $ 为速度张量； $ {x}_{i} $ 和 $ {x}_{j} $ 分别为空间上的分量； $ {\varOmega}_{ij} $ 和 $ {S}_{ij} $ 分别为速度梯度张量的对称和反对称部分。

当Q＞0时，表示流体微团的旋转克服了剪切应力并在流动中占主导作用，可以判断流场中产生了漩涡。涡量通常用来描述涡旋强度，定义为流体速度矢量的旋度。

图11为Q值与涡量的比较。从图中可以看到，瞬时涡量图和Q值图都显示了在多个接触面互相影响下产生了众多大小不一的涡，在不同接触面之间互相影响的区域使用Q准则能够判断识别出更多的涡，而使用涡量判断则效果不明显；涡量图中顶部流经中心接触面、中心接触面流经底部的区域显示有大量涡存在，而该区域实际上是氦气沿着球面流动，流动分离并不剧烈，只会出现少量旋涡，可见涡量判据出现了误判。误判的根本原因还是在于涡量对于旋转与剪切两种流动一视同仁，在多曲面问题上无法准确判断旋涡，这也是涡量的局限性所在。因此，在后续分析中本文将采用Q准则判断涡区域。

图12、13为中心燃料球表面速度场与温度场的时均结果与瞬时结果。在燃料球顶部和底部可见两块低速区域，该低速区域的形成是因接触面1的存在挤压了流道，阻碍了来流氦气的流动。所有接触面后沿都形成局部高温热斑，而接触面前沿则存在低温区域，这是因为接触面后沿流道较为狭小，氦气流经此处无法充分换热，而前沿处于流道较宽阔地带，热量能够被来流氦气带走。顶部接触面3与底部接触面4之间的垂直区域是低速高温区且变化不大，而接触面2与接触面3之间的三角区域、接触面2与接触面4之间的三角区域有较大的温度和速度变化。因此，基本可以判断中心接触面2的存在直接改变了氦气流向，从而影响了对流换热。

图14为燃料球外0.6 mm处时均Q值分布。由图中可以看到，涡强度较高的区域集中在接触面2和接触面3之间的三角区域，而燃料球顶部接触面2与底部接触面4处产生的涡强度较小。结合图12时均速度和时均温度分布可以发现，在高速区域涡的强度很大，而接触面后沿几乎无涡，速度几乎为零，低温区域则完全涵盖了存在涡的区域。由此说明：低速区域基本不会产生涡，高速区域必然会产生大量的涡；涡的存在有利于氦气的流动，高速氦气带走了热量从而降低了球面温度。

从能量方程对流项可知，除了速度大小，速度矢量与温度梯度夹角也决定了速度与温度梯度点乘的大小，且夹角越小对流换热越好。图15为燃料球外0.6 mm处时均速度与温度梯度夹角分布。结合图14、15可以发现，在顶部与底部涡分布较少的区域，时均夹角普遍大于90°，部分区域甚至达到120°，过大的时均夹角对换热不利；而有漩涡产生的区域时均夹角普遍小于等于45°，此时漩涡加强了对流，增强了换热。在接触面2与接触面3、接触面2与接触面4共同作用的区域，大量涡的存在使得时均夹角分布不均，此处漩涡不如顶部与底部漩涡的作用明显。

2.4 BCP结构与BCC、FCC结构对比

BCP结构与常见的BCC结构都为体心结构，但是前者接触位置更多，结构更为紧凑；BCP结构与FCC结构相比，两者均具有相同的空隙率，但几何构成却大不相同。参考文献［9］、［10］，本文对3种结构做一简单对比。

图16为3种规则排列球阵在相同工况、相同位置温度场的差异。观察轴向截面结果，BCC结构体心球与球之间距离近，从而导致球间流域温度高于BCP和FCC结构；FCC结构由于其面心排布，面心球之间距离和BCP结构体心球之间距离相近，因而在球间流域有相似的结果，氦气流动换热效果较好；BCP结构顶部燃料球之间互相接触产生的局部高温热点使整体温度分布更不均匀，而这种接触在另外两种结构中并不存在，这也是BCP结构最大的特点。再来观察径向截面结果，BCC结构中心球与其他中心球之间并无接触，间隙处存在4块低温区域；而FCC结构与BCP结构则存在接触，并产生了高温热点，热点后方则有8小块低温区域，这些低温区域代表氦气沿球面流动而未发生流动分离。接触面的存在打断了FCC和BCP结构球面的连续性，故而使流动更加复杂。

3 结　论

本文采用大涡模拟方法，对比了面接触与人工间隙处理方式的差异，对面接触BCP结构的速度场与温度场进行求解，得到流场中的瞬时和时均速度和温度分布，并了解了其流动换热的特点。同时充分发挥LES方法的优越性，对流场中的涡进行了分析。最后得到的主要结论为：

(1) 人工间隙处理方式对于高温区域存在误判，采用面接触方式更为合理。在接触面后沿出现显著高温热点，且随着来流方向接触面越多，温度越高。在燃料球顶部与底部滞止区与尾流区氦气分离再附着造成局部高温区域，且滞止区高温面积小于尾流区。

(2) BCP排列球阵四个不同位置的接触面构成的三角形区域存在大量附着与分离流动共存的现象。

(3) BCP结构顶部球之间的互相接触使得流域被多种接触面分割，与BCC和FCC结构相比，其温度分布更为复杂。

(4) 涡在接触面附近、接触面三角区域、接触面垂直区域会影响温度分布，且从整体来看涡的存在对流动换热是有利的。