M310堆顶冷却结构流场和温度场数值仿真
针对M310堆顶通风冷却系统，建立了包括控制棒驱动机构的堆顶冷却结构仿真模型，根据冷却结构气体流动特点，应用六面体网格进行区域离散，基于ANSYS CFX13.0软件，数值仿真了堆顶冷却结构内流场和温度场，分析了不同冷却空气流量、环境温度以及围板与磁轭线圈之间间隙对冷却系统性能以及速度不均匀系数的影响。计算结果表明，在给定的发热量下各磁轭线圈表面温度较低，最高不超过180℃；增大冷却流量，磁轭线圈表面温度降低，速度不均匀系数增大；增大围板与磁轭线圈之间的间隙，有利于降低各层线圈之间的速度不均匀系数。
北京航空航天大学 航空科学与工程学院 谢永奇 宋俊
反应堆堆顶结构是反应堆的重要部件之一，其重要功能是冷却控制棒驱动机构（Control rod drive mechanism， CRDM），保证CRDM运行时的温度满足设计要求。由于CRDM外形复杂，工作工况多变，对结构通风冷却性能影响的因素较多，开展堆顶冷却结构冷却性能的试验研究，不但耗时耗资，而且难以了解其详细信息。
近年来随着数值仿真技术的发展，CFD方法已经越来越多的用于系统尺度的反应堆热工水力问题的研究，而且其计算结果较好，接近试验结果。但是，对于堆顶结构通风冷却系统的研究相对较少，现有研究往往在大量简化的基础上进行，例如在系统尺度计算中将磁轭线圈简化，仅考虑阻力损失系数的方法。而对于CRDM磁轭线圈进行详细考虑的系统尺度研究的公开文献还未见有报道。
本文针对M310堆顶通风冷却系统，建立了考虑CRDM的堆顶冷却结构仿真模型，基于ANSYS CFX13.0软件，分析了不同冷却流量、环境温度以及结构尺寸等因素下的堆顶冷却结构流场和温度场，研究结果对于详细了解堆顶冷却结构的流动特点以及结构的优化设计具有重要意义和工程价值。
一、计算模型与方法
1.物理模型
M310堆顶结构通风冷却系统由抽风风机、风管、通风罩、风管过渡段、围板、仪表及电气部件等组成，用于为CRDM磁轭线圈提供满足设计要求的冷却风量，其结构如图1所示。
M310堆顶冷却结构采用上进风，冷却空气通过围板上部开口进入到CRDM磁轭线圈区域，在下风罩内气流分别进入沿周向布置的四个风管，然后经风管过渡段进入排风管道由抽风风机排出。
建模中，为尽可能减小模型简化对计算结果造成影响，不考虑小的凸台、圆角以及可能存在的缝隙，建立65根CRDM管及磁轭线圈壳体。为了减小计算误差，在围板上部扩展计算区域，在AD连接管和BC连接管连接处增设出口风管，所建立的计算域模型如图2所示。
在原方案基础上，提出了两种改进方案，即围板与磁轭线圈之间的间隙缩小为6mm（简称Opt1）和增大为20mm（简称Opt2），以考查结构变化对速度不均匀系数与冷却效果的影响。
2.数学模型
本文假设堆顶通风冷却结构中空气不可压缩，且为牛顿粘性流体，流动为稳态紊流，考虑压力容器对下风罩区域的漏热以及CRDM及磁轭线圈之间辐射换热，则堆顶冷却结构内定常不可压缩空气流动和传热的控制方程可描述如下。
连续性方程：
（1）
动量方程：
（2）
式中， 为速度矢量； μeff为有效粘性系数； P为修正压力。
能量方程：
（3）
式中，λ流体导热系数；ST为源项。
湍流动能k方程：
（4）
湍流动能ω方程：
（5）
式中， PK为湍流生成率， μt = ，常数的取值如下：σ=2， σω=2，C1=0.09，C2=5/9，C3=0.075。
3.网格划分与敏感性分析
M310堆顶冷却结构复杂，系统最大结构尺寸与最小结构尺寸比很大，这给网格生成和计算成本方面造成了许多困难。采用四面体单元虽然可以方便网格划分，但是单元质量难以控制，为获得满意的计算结果所需要的网格数量较多。而六面体单元比四面体单元具有更好的流动信息捕捉能力，其形状较为规则，网格节点数相对较少。考虑到冷却结构内部实际流动特性，合理安排网格的疏密分布，在磁轭线圈之间、下风罩区域、风管转弯等局部区域设置足够密的网格分布，以适应速度梯度过大的影响，因此，采用非均匀的六面体网格进行区域离散，如图3所示，这样既能完整描述流场特征，又能减少划分网格数量，降低对计算资源的需求，提高计算速度。
由于磁轭线圈区域（上风罩）和下风罩区域内气流流动是考察重点，这两部分网格数量划分较多。对如表1所示的4种网格划分方案进行网格敏感性分析。需要说明的是，上风罩区域CRDM磁轭线圈之间的网格单元层数和沿CRDM长度方向的分段数对计算结果会产生较大影响。由于上风罩围板和磁轭线圈壁面之间间距较小，应当给出不少于5层的网格单元，并考虑第一层单元距离壁面的高度大小。沿CRDM长度方向上，网格节点数不少于121个。
针对表1给出的不同网格方案，计算得到的上风罩区域进出口的压差分别为-416.2Pa、-422.8Pa、-422.6Pa和-424.3Pa，可以看出在上风罩区域网格数量达到532万时，上风罩进出口压差值与网格数量为671万时的差别很小，因此，当区域划分加密得到的网格单元数大于1500万个时（其中磁轭线圈区域网格量需大于530万），计算结果相对于该单元数情况下的误差小于2%，故可认为此时计算结果已与单元数的多少无关。各方案节点数和体单元数相差不大，节点数1700万个， 体单元数1500万个。
4.边界条件与求解策略
计算中为更好反映冷却空气温度与密度之间的对应关系，采用理想气体模型作为实际空气的近似。堆顶通风冷却系统流动和传热边界主要包括气流入口、气流出口和各个壁面。假定环境压力为101325Pa，环境温度即为进气温度，进口边界条件设置为总压，压力为101325Pa，出口边界条件设置为质量流量。冷却结构各壁面的流动边界条件采用固体壁面无滑移边界条件，应用壁面函数决定靠近壁面附近流体的流动行为。根据所给定的每根CRDM的发热量为10kW，CRDM及磁轭线圈壁面设定为热流密度边界条件。反应堆压力容器通过保温层向下风罩内漏热，设定保温层壁面为恒壁温边界，温度为70℃。不考虑其它壁面的传热，即设置为绝热壁面。
采用标准k-ω模型封闭湍流控制方程，Thermal Energy模型模拟传热，Discrete Transfer模型模拟辐射换热，并假定各个方向散射相同。控制微分方程离散时采用有限体积法，应用High Resolution方法求解连续方程、动量方程、能量方程及湍流方程，时间尺度控制选择物理时间尺度，各方案计算至200步左右达到收敛。
二、仿真结果与分析
1.原方案结果分析
图4为原方案冷却结构速度分布和局部放大图。可以看出，冷却空气在上风罩进口出的气流流速较低，大部分不超过5m/s；冷却空气进入CRDM磁轭线圈间隙后，在风机负压抽吸作用下，气流流速增大，大部分不低于14m/s。气流经围板区域进入下风罩区域后，由于流动空间增大，气流流速有所降低；然后经四个连接管弯头后进入AD连接管和BC连接管，管内的气流流速较高，大部分不低于25m/s，在局部区域（吊杆穿过连接管处）气流流速超过50m/s；两股气流汇合后进入出口风管。由图4可知，系统内存在较多的漩涡，这将导致较大的压力损失。
图5分别给出了原方案z=0m的 XY平面上位置示意及其上气流速度、压力分布云图。由图5（b）可以看出，在CRDM之间间隙内，气流流速较高（中间红色区域），大部分不低于18m/s，在上下磁轭线圈之间，气流流速也存在小流速区，大部分不超过3.0m/s。两侧流道内的气流流速要低于中间区域的流速。气流流经线圈后，气流流速因流道截面变大而减小，同时受下风罩内负压影响，下风罩内气流流动方向偏向连接管与下风罩接口处。
由图5（c）可以看出，在下风罩区域内静压大部分在-370Pa～-480Pa，沿CRDM杆向下，流动内静压逐渐减低，在靠近最下层线圈局部区域，由于漩涡的存在其压力也较低，该截面上的静压最低为-532.9Pa。
图6分别给出了上、中、下三层线圈中间截面上的压力分布，对应截面为y=1.58m、 y=1.25m、 y=0.95m的 ZX平面。可以看出，最上层的静压最大，中间次之，最下层静压最小，最上层和中间层静压分布较为均匀，最下层静压分布差别较大，中间区域的静压在-444Pa左右，四周的静压在-530Pa左右。总的来看，上层和中层线圈间区域靠近围板壁面处气流静压较大，但在下层靠近上风罩中间区域静压较大，这是由于风道负压造成的。
图7给出了连接管汇合处截面上的速度分布，可以看出，在出风口汇合处，来自AD和BC连接管的两股气流以较高的流速掺混，在该截面上相互干扰，形成中间流速小、两侧流速高的速度分布，两侧流速大部分不低于30m/s，中间流速相对较低，大部分不超过21m/s。结合图4可知，该位置存在漩涡。这也说明了在出口扩展计算区域是必要的，如果直接以该截面作为出口截面，施加压力边界条件，将会带来较大的计算误差。
无论冷却系统是送风还是吸风，流场不均匀系数有两种计算公式：
其中，v大、v小、 分别为间隙中心测得的各点风速的最大值、最小值和平均值，一般的有ε大<ε小，所以应采用ε小作为评定流场均匀性的标志性参数。
表2给出了原方案上、中、下三层线圈每层中间截面上的速度不均匀系数。由表2可知，原方案上层线圈的不均匀系数为19.2%，中间线圈的不均匀系数为17.5%，下层线圈的不均匀系数为22.7%，均大于所要求的20%。
图8给出了原方案CRDM的温度分布云图，由于上风罩进口的气流速度较低，在上风罩入口区域的CRDM杆上的温度较高，但是上、中、下三层线圈壁面的温度大部分不超过130℃，个别区域的温度不超过213℃，最上层线圈壁面温度较高，中、下层壁面温度较低，其温度分布与相应的流场分布对应。根据流场分析结果可知，相邻磁轭线圈之间的气流流速均较高，且从上到下气流流速有所增加，所以中、下层线圈的壁面温度较低。由于在上下层线圈之间气流存在漩涡，且流速不高，所以这部分区域壁面温度相对较高。
2.不同冷却风量结果分析
在原方案计算中风机体积流量为70600m 3/h，在此基础上通过增大流量5%和20%，以考查冷却结构的流场分布特性。
当增加的风机风量为原来20%时，整个冷却结构的冷却空气流量为26.2kg/s，出口空气的温度为34℃。风量增大5%和20%时冷却空气沿流道流动过程中的压差如表3所示。由表3可以看出，风量增大，围板区域的流动阻力增大，进出口压差为-648.5Pa；AD和BC连接管内的气流流动阻力损失为-1 209.7Pa和-1 240.4Pa，两个管道内压力差基本相同，对应的空气流量分配也相同。整个模型进出口的压力降为2996Pa，相对于流量为70 600m 3/h的压降均增大。
对于风量增大5%，图9给出了z=0m的XY平面上的速度矢量分布图（包括进气口附加区域、围板和下风罩区域）。由图9可以看出，在CRDM之间间隙内，气流流速较高，大部分不低于18.2m/s，在上下磁轭线圈之间，气流流速也存在小流速区，大部分不超过3.0m/s。两侧流道内的气流流速要低于中间区域的流速。气流流经线圈后，气流流速因流道截面变大而减小，同时受下风罩内负压影响，下风罩内气流流动方向偏向连接管与下风罩接口处。在下风罩区域内，气流应存在多个漩涡。
在靠近围板上端（图9中局部放大图），外部气流受围板壁阻挡，绕流进入上风罩内，并在靠近围板壁区域形成漩涡，这会对靠近围板壁间隙内的气流流动造成不利影响。
增大流量20%后，上、中、下三层线圈平面上的速度不均匀系数分别为20.5%、15.8%和25.3%。相对于未增大时均增大。而流量增大5%时，各层线圈之间的速度不均匀系数略有增加。
风量增大20%时整个CRDM线圈和杆的温度分布与原方案的相似，CRDM壁面的最高温度为182.5℃，高温区主要位于上风罩入口区域的CRDM杆上。但是上、中、下三层线圈壁面的温度大部分不超过101℃，最上层线圈壁面温度较高，中、下层壁面温度较低，其温度分布与相应的流场分布对应。与原方案相比，CRDM温度降低了约30℃。
3.不同环境温度结果分析
该反应堆所在地域的环境温度为5℃～55℃。环境温度的变化将导致空气密度的变化，进而影响风机质量流量。计算中其他条件设置与原方案相同，在环境温度为5℃时，冷却系统空气流量为22.9kg/s，出口空气温度为23.4℃。在环境温度为55℃时，冷却系统空气流量为19.6kg/s，出口空气温度为73℃。
图10给出了55℃时下风罩区域 y=0.5m的 ZX平面上速度和压力分布，可以看出，在下风罩中心区域（5个CRDM），气流流速较低，四个通风弯头之间区域的气流流速也不高，大部分不超过4.0m/s。在靠近通风弯头区域的气流流速较高，不低于11.6m/s，弯头内部分气流流速超过了18m/s。从速度分布上看，A、B通风弯头内气流速度（平均速度为34.4m/s和34.1m/s）高于C、D通风弯头内气流速度（32.1m/s和31.1m/s），这也可能引起上风罩区域速度分布的不均匀。由图8（b）可以看出，下风罩中间区域静压分布在-364.0Pa～-431.7Pa，平面静压不均匀值超过60Pa。靠近通风弯头的静压较低，静压大部分在-517Pa～-635Pa。
在环境温度为5℃时，冷却风量变化较小，各层速度不均匀系数变化不大。在55℃时，冷却风量相对于原方案减小了15%，相应的上、中、下三层线圈平面上的速度不均匀系数分别为18.2%、13.3%和20.5%，相对于原方案均有所改善。但是，由于冷却空气进口温度升高，相应的CRDM线圈的温度也升高，局部温度超过了165℃。
4.不同改进方案结果分析
图11分别给出了Opt1和Opt2方案 z=0m的 XY平面上位置示意及其上的速度矢量分布图（包括进气口附加区域、围板和下风罩区域）。由图11（a）可以看出，在CRDM之间间隙内，气流流速较高，大部分不低于16.1m/s，在上下磁轭线圈之间，气流流速也存在小流速区，大部分不超过3.0m/s，而且存在两个对称的漩涡。两侧流道内的气流流速要低于中间区域的流速。气流流经线圈后，气流流速因流道截面变大而减小，同时受下风罩内负压影响，下风罩内气流流动方向偏向连接管与下风罩接口处。在下风罩区域内，气流应存在多个漩涡。
由图11（b）可以看出，在CRDM之间间隙内，气流流速较高，大部分不低于15.6m/s，调整间隙值后，CRDM之间气流流速最大值由Opt1的29.0m/s降低为23.3m/s，速度分布与Opt1方案相似。在进气口附加区域，流线较疏，对应的气流流速较小。外界环境空气绕流围板壁进入上风罩区域，在靠近围板与线圈之间间隙区域，气流出现了漩涡，如图11（b）中局部放大图，这不利于靠近围板壁间隙内的气流流动。但相对于Opt1的改进方案，该处漩涡区域有所减小，而且漩涡靠近围板上端。
图12分别给出了原方案、Opt6方案和Opt20方案上风罩进口截面上的速度分布，可以看出，Opt6、Opt20方案进口截面上的速度分布与原方案相似，也是在靠近围板壁面处速度相对较高，中间区域速度分布均匀。所不同的是，减小围板与线圈之间间隙，使得靠近围板处的速度最大值略小于原方案。增大围板与线圈之间间隙，有助于增大间隙内气流流速，使得截面上的速度分布更趋于均匀。
表4和表5分别给出了两种改进方案的速度不均匀系数值。由表4可知，上、中、下三层磁轭线圈其速度不均匀系数分别为24.0%、24.4%和18.0%。减小围板与线圈之间的间隙，上层和中层线圈间的速度不均匀系数均增大，下层线圈间的速度不均匀系数减小，这可能是由于间隙的减小增大了下层线圈周围的气流流速所致。这表明减小围板与线圈之间的间距，不利于改善CRDM间的流场不均匀性。
由表5可知，上、中、下三层磁轭线圈其速度不均匀系数分别为16.5%、6.9%和18.8%，满足设计要求，这表明增大围板与线圈之间的间隙，有利于改进流场不均匀性。
图13给出了改进方案Opt1的CRDM温度分布，可以看出，CRDM壁面的最高温度为228.3℃，高温区主要位于上风罩入口区域的CRDM杆上，上、中、下三层线圈壁面的温度大部分不超过135℃，最上层线圈壁面温度较高，中、下层壁面温度较低，其温度分布与相应的流场分布对应。
改进方案Opt2的CRDM温度分布于Opt1相似，其最高温度为208℃。综合分析两种改进方案与原方案可知，增大围板与线圈之间的间隙，上、中、下三层线圈上的速度不均匀系数均减小，尤其以中间层线圈不均匀系数减小最多，CRDM磁轭线圈壁面温度也降低，表明增大围板与线圈之间的间隙有利于改善其流场均匀性和通风冷却效果。
三、结语
本文应用ANSYS CFX13.0软件针对M310堆顶冷却结构流场和温度场进行了数值仿真研究，主要得到如下结论。
（1）增加冷却空气流量，系统流动阻力增大，CRDM温度降低，CRDM间隙中气流流速增大，上、中、下三层线圈间的气流不均匀性也增大，其中下层线圈受影响较大。
（2）环境温度升高，冷却空气质量流量减小，上、中、下三层线圈间的气流不均匀性变化不大，CRDM壁面温度增大，不利于系统冷却。
（3）减小围板与线圈之间的间隙，不利于降低线圈间速度不均匀系数，增大围板与线圈之间的间隙，有利于提高线圈间速度均匀性，增强CRDM的冷却效果。
本项目研究结果可为M310堆顶冷却结构的优化设计及工程实施提供科学依据。