李小畅，邓 坚，杨小磊，郭涌辉，陈广亮，田瑞峰,*

(1.哈尔滨工程大学 航天与建筑工程学院，黑龙江 哈尔滨 150001；2.中国核动力研究设计院 核反应堆系统设计技术重点实验室，四川 成都 610041；3.哈尔滨工程大学 核科学与技术学院，黑龙江 哈尔滨 150001)

反应堆三维精细化数值建模与仿真技术是克服传统研究方法各种不足及大幅降低研发成本和周期的有效途径之一，也是实现核能技术多学科耦合协同研发与设计的重要手段，因此受到世界各国重视。各核电大国已投入大量资源进行相关研究并组建了相应核能数值仿真研究中心，并确定了具体研究项目[1-2]，如美国的CASL和NEAMS项目，欧洲地区的NURESIM系列项目等，中国核动力研究设计院联合哈尔滨工程大学参考美国CASL项目确定了“数值反应堆”研发项目。国内也加大了推进核能一体化数值建模与仿真技术开发力度，以便在新一轮核能研发领域赶上欧美发达国家的步伐。反应堆一体化数值建模与仿真技术研究对象往往是设备级的，如全尺寸或1/2尺寸反应堆，且需考虑多物理场耦合求解，如热-流-固-物理耦合。反应堆三维精细化多物理场耦合仿真的前提是解决各单一物理场的三维全尺寸级精细化求解难题，当前仍有诸多基本科学与技术问题尚待明确和解决，本文所研究的反应堆燃料组件三维精细化热工水力数值仿真求解效率问题则是其中之一。

压水堆堆芯包含数十至两百多个燃料组件不等，典型燃料组件由17×17的燃料棒方形排列而成，同时在竖直方向上设置约8个起固定和搅混作用的搅混翼定位格架，总高度约4 m。初步估算，即使基于计算量较小的雷诺时均+湍流涡黏模型+多面体网格的计算方法，对单个完整燃料组件进行三维精细化热工水力数值求解，所需计算单元总数约为1亿；采用大涡模拟(LES)结合壁面函数法则需至少10亿网格单元[3]。尽管近20年来，国内外针对CFD方法应用于反应堆燃料组件的相关研究成果层出不穷，且已证明是一种成熟有效的方法[4-8]，然而这一阶段三维精细化热工水力分析大多存在一个共同特点：仿真对象主要基于简化模型或局部模型，如3×3、5×5等，与典型真实燃料组件17×17的规模还存在较大差距，主要原因一方面是早期研究处于验证分析阶段，采用小规模棒束可提高研究效率；另一方面则是受到计算能力与计算效率的限制。已有研究表明，真实燃料组件的热工水力性能无法完全用小规模棒束通道来表征，不同棒束规模燃料组件热工水力特性存在明显差异[9-10]。可见，目前无论是针对新一轮核能技术研发、核安全分析及由此衍生的“虚拟反应堆”研发，还是针对反应堆燃料组件热工水力特性理论研究本身，都对反应堆燃料组件大规模棒束通道的三维精细化数值仿真有迫切需求。然而当前面临的一个非常现实的问题是，如何解决反应堆全尺寸级三维精细化数值仿真所面临的网格数量过多、数值计算量过大及积累误差过大等问题。

对于压水堆燃料组件热工水力数值求解，造成其数值计算量过大的根本原因在于基于欧拉法与贴体网格技术的三维精细化数值分析方法不得不采用数量极为庞大的“非结构化网格”来表征复杂结构的真实物理边界，导致计算成本过高。以往针对这一问题的主要解决方法为采用多孔介质模型或粗网格技术等降低计算量[11-14]。直到近年来，学术界才陆续出现专门针对反应堆堆芯的三维精细化数值仿真计算量与求解效率的相关研究。Hu等[15-16]首先对钠冷快堆中的绕丝棒束通道进行相关研究，提出了一个动量源数理模型来表征绕丝对棒束通道热工水力的影响。Capone等[17-18]针对压水堆燃料组件，提出了一个基于高分辨率湍流模型(如大涡模拟)计算结果的动量源模型，并将此方法应用到一5×5的棒束通道中。Mikuz等[19]同样基于燃料组件完整贴体网格模型的计算结果，提取了相应的动量源项，并在一3×3的裸棒束通道中实现了搅混翼格架对流场与温度场影响的模拟。上述针对压水堆的相关研究，大多基于完整燃料组件的贴体网格计算结果来提取并拟合出合适的动量源模型，本文则直接基于复杂结构的影响机理建立其数学模型。

本文主要基于燃料组件搅混翼对流场及温度场的影响机理，将复杂结构对热工水力特性的综合影响归结为流体质点受到的某种虚拟体积力，只需给出这种力的数学描述并以动量源项的形式加入动量方程，即可在无翼片通道中获得搅混通道的热工水力特性。通过与实验数据及传统贴体网格建模计算结果的对比分析，验证本文所建立虚拟体积力动量源模型的有效性。

1 虚拟体积力动量源建模

本文以Karoutas等[4]的搅混翼棒束通道激光多普勒测速(LDV)实验为研究对象，建立压水堆燃料组件搅混翼的虚拟体积力动量源模型。该实验提供了1个5×5棒束通道搅混翼格架下游不同位置的流速测量结果，实验模型及尺寸如图1所示。真实压水堆燃料组件的搅混翼格架，除搅混翼外，还包含对燃料棒起约束作用的刚凸及弹簧，本文所研究的虚拟体积力动量源模型主要针对搅混翼，暂不考虑刚凸及弹簧。

图1 5×5棒束实验模型及两通道计算模型Fig.1 5×5 rod bundle experiment model and two-subchannel model

典型燃料组件的搅混翼具有周期性重复排列的特点，最小周期性单元为图1中蓝色虚线框所示的四通道模型。对于燃料组件的中心区域棒束通道，已有研究结果[10]表明，在数值建模与仿真中，四通道模型最小周期性单元还可采用图2所示的周期性边界交叉配对(P1a-P1b、P2a-P2b、P3a-P3b)的方式进一步简化成两通道模型而不丢失其主要热工水力特征，进一步提高研究过程中的数值求解效率。对于带搅混翼的棒束通道模型，由于每个狭长的子通道内至少存在2个偏折方向不同的搅混翼，因此在对搅混翼棒束通道进行数值求解时，基于贴体网格理论的CFD方法需将包含搅混翼格架的流体域单独隔离出来进行非结构化网格建模，其他规则区域则将面网格拉伸成体网格以减少网格数量和提高网格质量，如图3所示。可见，搅混翼格架的存在是造成贴体网格建模数值计算量过大的根本原因。反应堆堆芯燃料棒的排列本身非常规则，其数值分析可采用高质量结构化网格建模。如果能通过某种数学物理模型来描述搅混翼对棒束通道内流场与温度场的影响，则可将复杂搅混翼棒束通道的数值建模简化为无翼片通道建模，从而大幅降低数值计算量。本文以图1中的两通道模型为研究对象，建立搅混翼的虚拟体积力数学模型，以考虑搅混翼对流场与温度场的影响。

图2 两通道模型周期性边界Fig.2 Periodic boundary conditions of two-subchannel model

图3 燃料组件搅混翼棒束通道典型几何与网格模型Fig.3 Typical geometrical and mesh models of fuel assembly rod bundle with mixing vane

通过理论分析可知，搅混翼对棒束通道流场与温度场的影响主要体现在：搅混翼所产生的二次流及旋涡增强了流体的横向运动及不同子通道间的交混能力，从而强化了流体与壁面间的对流换热并降低了温度不均匀性，同时增大了压降。进一步分析可知，换热的强化及压降的增大本质上均由搅混翼所产生的二次流及漩涡的交混效应所造成，而二次流与漩涡运动又主要由搅混翼的阻流和导流作用引起，棒束通道下游典型横向流场及涡结构如图4所示。因此，可认为导致上述各物理现象的原始驱动力是由于搅混翼对流场中的流体质点施加了某种作用力，这种作用力改变了流体的运动规律并导致了复杂的流动现象(二次流、漩涡、搅混等)及换热现象(降低温度不均匀性、强化换热)。因此，如果能基于搅混翼对流场与温度场的作用机理给出这种力的数学描述，则在真实燃料组件的数值建模中就无需再对搅混翼进行数量庞大的非结构化贴体网格建模，从而简化建模过程及大幅降低计算量、提高计算效率，最终实现在简单规则的裸棒束模型通道中即可获得搅混翼棒束通道中的复杂流场，而其他所有由流场变化引起的物理现象，如换热的强化、温度不均匀性的降低、压降的增大等均可自然获得。

图4 搅混翼格架下游横向流场及流型示意图Fig.4 Lateral streamline and flow pattern downstream of mixing vane grid

搅混翼虚拟体积力示意图如图5所示。假设棒束通道中任一搅混翼对流体的虚拟体积力为Fb，结合翼片形状可将其分解为3个方向的分力，即与搅混翼垂直的法向分力Fn、与搅混翼平行的分力Fp、与搅混翼相切的分力Ft，则该搅混翼的虚拟体积力合力如式(1)所示。还可进一步将搅混翼各方向的虚拟体积力在直角坐标系下分解，以便后续数值求解中动量源项的建立及加载，如式(2)所示。

Fb=Fn+Fp+Ft

(1)

(2)

式中，θ为搅混翼相对燃料棒轴向的偏折角，压水堆燃料组件约为25°。结合实际搅混翼对流体的阻流、导流及搅混现象，从图5可看出，垂直于搅混翼平面的力Fn一方面对流体流动产生阻碍作用，同时在x正方向提供一定程度的横向驱动力；平行于搅混翼方向的力Fp主要提供额外的摩擦阻力，同时对流体在x轴负方向提供一定驱动力；与搅混翼相切的力Ft则主要影响流体在y方向的流动。这些力既能把流体导向某个方向并加速，也能对流体产生一定的阻碍作用，其综合效果体现了搅混翼对流场的影响。

图5 搅混翼虚拟体积力示意图Fig.5 Schematic of virtual body force for mixing vane

根据式(3)所示的不可压缩黏性流体动量守恒方程，搅混翼的虚拟体积力将以动量源项的形式加入到动量方程的体积力项f(N/m3)中，其物理意义为单位体积流体所受到的体积力。

(3)

式中，ρ、V、p、t分别为密度、流速、压强及时间。

在搅混翼法向方向，流体质点受到搅混翼的阻碍及导向作用后的加速度可由流速对时间的质点导数(全导数)来表示：

(4)

式中，nn、nt、np分别表示搅混翼的法向、平行及切向方向。稳态工况下，式(4)右侧的当地加速度项∂Vn/∂t=0，只剩迁移加速度项。考虑翼片法向的体积力Fn是引起流体在该方向上动量变化的主要原因，则法向虚拟体积力导致的控制体单位体积内流体动量的变化率可表示为：

(5)

式中：Vn、Vt、Vp分别表示翼片法向、平行及切向方向的流速；nn为法向单位矢量。对于式(5)中法向流速在各方向的速度梯度，本文结合搅混翼几何特征，采用式(6)所示形式。

(6)

式中：k为乘数因子，可基于具体燃料组件的模型实验数据对计算结果进行修正；ln为搅混翼法向厚度；at、ap分别为切向厚度修正因子及平行厚度修正因子。考虑到搅混翼实际导流效果，法向速度在翼片平行方向的梯度必定大于在翼片切向方向的梯度，因此需保证at>ap。

对于与搅混翼平面平行和相切的两个方向的虚拟体积力Fp及Ft，考虑到这两个方向的力主要由翼片与流体之间的切向力引起，因此基于达西公式建立其动量源模型：

(7)

式中：np、nt分别为搅混翼平行向及切向单位矢量；ζp、ζt分别为各自方向上的阻力系数，由McAdams关系式[20]确定。

ζ=0.184/Re0.23×104

(8)

由式(5)～(8)可看出，搅混翼对流体各方向的虚拟体积力均可表示成当地速度的函数。数值求解过程如图6所示，首先基于无搅混翼的棒束通道数值模型，根据初始流场计算原搅混翼所在空间网格节点各方向的虚拟体积力动量源项。将所计算的动量源项加载到动量守恒方程中并迭代，获得新的流场后判断结果是否收敛，若未收敛，则基于新的速度场重新计算动量源项，重复上述过程直至计算收敛。虚拟体积力动量源项与当地速度场之间为自动耦合关系，基于当地速度得到动量源项，动量源项再作用于流场，直到所计算的动量源项恰好能维持该速度场的稳定。可见，整个动量源计算及求解过程只需基于当地速度，无需引入其他变量。虚拟体积力模型在获得搅混翼下游流场的同时，翼片对温度场的影响可自然获得。这是因为，所形成的横向流场将一个通道的冷却剂输送至另一个通道降低了温度分布不均匀性；各通道内的漩涡增加了流体的搅混效应，并迫使边界层内的高温流体向主流低温区运动，提高了加热壁面与流体间的温差，实现了换热的强化。

图6 数值求解过程Fig.6 Numerical solution process

图7 搅混翼类型Fig.7 Type of mixing vane

从上述建模过程可看出，各方向的虚拟体积力动量源数学模型的具体形式及其在直角坐标系下的分解与翼片位置及偏折方向有关。整个反应堆堆芯燃料组件的搅混翼是4种类型翼片(图7)周期性重复排列的结果，与坐标轴相对位置及偏折方向相同的搅混翼具有相同形式的虚拟体积力动量源项。因此，只需对4种类型的搅混翼分别建立其虚拟体积力动量源模型，并加载到同类型的翼片上，即可实现完整燃料组件的数值建模。另外3种类型的搅混翼虚拟体积力动量源建模方法与翼片类型1类似。

2 数值验证

基于Karoutas等[4]的搅混翼棒束通道激光多普勒测速(LDV)实验，以图1中两通道模型为研究对象，分别建立其虚拟体积力动量源数值模型及贴体网格数值模型，其中侧边界采用图2所示的交叉配对周期性边界条件，以考虑该计算通道与周围通道之间的交混。该实验来流速度为6.79 m/s，温度为26.67 ℃，系统压力为48.3 kPa。为研究虚拟体积力动量源模型对温度场的模拟效果，计算中对燃料棒施加707 kW/m2的定热流密度。数值求解基于ANSYS FLUENTL，分别建立两通道模型中4个搅混翼的虚拟体积力动量源数学模型，通过用户自定义函数(UDF)以动量源项的形式加入动量方程中。通过对比分析虚拟体积力动量源模型计算结果、常规贴体网格计算结果及实验数据来验证模型的有效性。

对于搅混翼格架棒束通道的网格建模，当搅混翼部分的网格基准尺寸约为0.25 mm时，可避免网格尺寸对计算结果的影响[10]，当采用SSTk-ω湍流模型结合y+为30左右的边界层网格时，对流场、压降及温度场有较好的计算效果[8]。基于两通道模型，分别建立带搅混翼的完整贴体网格模型及无搅混翼的混合网格模型，前者用于传统贴体网格求解，后者用于虚拟体积力动量源求解。从虚拟体积力动量源数学建模过程可知，所建立的动量源模型须加载到原搅混翼所在区域。鉴于本文主要是验证虚拟体积力动量源数学模型的有效性，在验证计算时首先将无搅混翼模型的流体域分为两部分，一部分为搅混翼所在区域的流体网格，另一部分为其他区域流体网格，以方便动量源模型加载到原翼片所在区域。经虚拟体积力动量源计算结果与实验数据及贴体网格建模结果的综合对比分析后，确定虚拟体积力动量源模型参数如下：k=2.2，at=1.3，ap=1.7。

采用传统贴体网格建模及虚拟体积力动量源建模所计算的横向流速与实验数据的对比示于图8，数据均取自图1中的实验测量路径。图中，Vy为y方向的流速，Vavg为取值路径上的平均流速，相对位置为x坐标值与棒间距P的比值，z为所在路径与搅混翼格架末端的轴向距离。从图8可看出，基于贴体网格建模与虚拟体积力动量源建模所计算的横向流速在分布趋势及大小上均与实验结果基本吻合，仅在搅混翼格架近下游区域的2个侧端面处与实验结果存在一定偏差。从动量源法与传统贴体网格法的对比结果看，在搅混翼格架的近下游区域，二者非常接近，远离格架则出现一定偏差，动量源法计算的横向流速衰减略偏大。

格架下游的横向流场本质上是搅混翼所引起的二次流现象，体现了搅混翼对下游流场的搅混与导流效果，也是换热强化、压降增大的根本原因。为更直观地展示动量源模型对搅混翼下游横向流场的模拟效果，采用两种方法计算了格架下游不同位置横向流场，结果示于图9。从图9可看出，两种方法所计算的搅混翼格架下游横向流场与图4所示的典型横向流型一致，通道中心存在1个大尺度漩涡，流体沿着燃料棒边缘从一个通道流向另一个通道，且在通道出口处存在小尺度漩涡。动量源法所计算的通道中心大尺度漩涡位置及形状与传统贴体网格计算结果基本一致，但二者在z=12.7 mm截面燃料棒近壁面区域的小尺度漩涡计算结果存在一定差别，贴体网格计算结果为单个小尺度漩涡，动量源法结果则分裂为2个，且位置略有偏差。图中各截面云图为轴向流速分布，可见动量源法计算的轴向流速的大小及分布与贴体网格结果也基本吻合。

图8 搅混翼格架下游横向流速分布Fig.8 Lateral velocity profile downstream of mixing vane grid

图9 搅混翼格架下游横向流场对比Fig.9 Comparison of lateral flow fields downstream of mixing vane grid

搅混翼格架及其下游压降计算结果示于图10，图中实验关联式结果基于Chun等[21]的半经验公式，该压降关联式由格架形状阻力、格架摩擦阻力、搅混翼形状阻力及燃料棒摩擦阻力4部分组成。从图10可看出，动量源法与传统贴体网格法所计算的各段压降与实验关联式均较为接近，局部最大偏差在10%以内。动量源法与贴体网格法计算的压降与实验关联式的偏差主要出现在搅混翼格架段，格架下游区域的摩擦阻力段则十分接近，但与实验关联式相比均略偏高。对比两种计算方法的结果可知，二者压降曲线在搅混翼格架下游部分重合度较高，说明动量源法很好地捕捉到了格架下游的阻力特性。

图10 搅混翼棒束通道压降Fig.10 Pressure loss along rod bundle with mixing vane grid

贴体网格法及虚拟体积力动量源法所计算的Nusselt数(Nu)沿燃料棒轴向及周向的分布规律示于图11。Nu的物理意义为壁面附近流体的无量纲温度梯度，其大小表征了流体与固体壁面间对流换热的强弱。从图11可看出，动量源法所计算的Nu沿轴向的分布与传统贴体网格计算结果吻合良好，最大偏差在5%以内。周向分布上，二者在搅混翼格架的近下游区域(12.7 mm)吻合较好，在远下游区域(50.4 mm)则存在一定偏差。远下游区域的偏差主要体现在Nu的局部极值在方位角上有所偏移，且该偏移与图8中的横向速度偏差存在一定关联性，原因在于Nu的局部大小主要受当地流速的影响。总体上，虚拟体积力动量源法所计算的Nu无论在数值还是分布规律上，均与传统贴体网格法计算结果吻合较好。从本文动量源数理建模及计算过程来看，在未对温度场求解进行任何特殊处理的情况下较好地捕捉到了搅混翼格架对棒束下游传热特性的影响，说明本文所给出的搅混翼对棒束通道流场与温度场的影响机理分析是合理的，也验证了虚拟体积力动量源建模思想的可行性。

为进一步综合评估虚拟体积力动量源模型对搅混翼棒束通道下游流动与传热特性的数值模拟效果，对搅混翼格架下游燃料棒近壁面对数律区流速沿轴向分布进行对比，结果示于图12。其中Dh为棒束通道水力直径。近壁面对数律区是边界层中的重要区域，该区域流速对固体壁面的对流换热起决定性作用，其大小与分布体现了搅混翼对下游流场与温度场的影响。从图12可看出，基于虚拟体积力动量源法所计算的燃料棒表面对数律区流速，无论是分布趋势还是大小均与完整贴体网格建模结果基本吻合。虚拟体积力动量源所计算的近壁面对数律区流速虽总体上略高于传统贴体网格计算结果，但最大偏差小于5%。此外，对比图11可知，Nu与对数律区流速的轴向分布存在一定相似性，符合流动与换热关系的理论预期。

图11 燃料棒轴向及周向Nu分布对比Fig.11 Comparison of axial and azimuthal variations in Nusselt number

为验证虚拟体积力动量源模型对不同来流工况的适应性，对不同流速(V0)工况下棒束通道内横向流速及对流换热系数与传统贴体网格建模法进行对比，结果示于图13、14。从图13可看出，不同来流工况下，两种方法所计算的格架下游12.7 mm处的横向流速均十分接近。可见从流场的角度，虚拟体积力动量源模型在各流速工况下均具有较好的适应性。从图14可看出，两种方法的计算结果在小流速工况下吻合更好，随着流速的增大，二者偏差呈现略增大的趋势。总体上，不同来流工况下虚拟体积力动量源法与传统贴体网格建模法所计算的对流换热系数在大小及分布趋势上均有较高的吻合度。

图12 燃料棒近壁面对数律区速度分布Fig.12 Profile of velocity in log-law region near rod wall

为明确虚拟体积力动量源法对棒束通道全结构化网格模型的敏感性，以评估计算效率的提升效果，本文进一步对两通道模型进行了六面体全结构化网格建模。求解时先标记原搅混翼所在空间的网格单元，用于动量源模型的加载。计算表明，当采用全六面体结构化网格建模，两通道模型网格单元数达到10.8万时，动量源计算结果不受网格影响，与贴体网格的对比如图15所示。

基于结构化和非结构化网格建模时动量源法所计算的当地横向速度及Nu轴向分布与带搅混翼贴体网格结果的对比示于图16。从图16可看出，当采用全结构化网格时，动量源法所计算的横向速度及Nu均与非结构化网格动量源模型及贴体网格模型整体吻合良好，均捕捉到了各参数的主要变化特征。仔细观察横向速度对比结果可知，结构化网格动量源模型所计算的搅混翼格架下游12.7 mm处的横向速度在左侧通道中较另外两种方法偏低，速度变化更趋平缓。可推测，该偏差产生的原因在于全结构化网格模型中施加动量源的网格单元无法与真实结构的边界完全贴合，丢失了部分结构特征对流场的影响。上述分析表明，基于全结构化网格的动量源法在较好地预测到搅混翼格架下游流动与传热特性的同时，可显著提高棒束通道的数值求解效率。以本文单跨两通道模型为例，采用传统贴体网格建模所需单元数为43.2万，基于结构化网格的动量源法则只需10.8万，网格数量减少了76%。可预测，将动量源法应用于单个完整燃料组件或堆芯，计算效率的提升将更可观。

图13 z=12.7 mm处不同流速工况下横向速度计算结果对比Fig.13 Comparison of calculated lateral velocities under different flow rates at z=12.7 mm

图14 不同流速工况下换热系数计算结果对比Fig.14 Comparison of calculated heat transfer coefficients under different flow rates

图15 用于动量源模型的结构化网格与用于完整结构的贴体网格Fig.15 Structured mesh for momentum source model and body-fitted mesh for full structure

图16 基于不同网格的动量源模型与贴体网格模型结果对比Fig.16 Comparison between results of momentun source model based on different meshes and body-fitting mesh model

3 结论

1) 压水堆燃料组件棒束通道中的搅混翼对流动与传热特性的综合影响，可由复杂结构对流体质点产生的某种虚拟体积力来表征，只需对这种体积力建立合适的数学模型并以动量源形式加入无搅混翼棒束通道模型中，即可同时获得搅混翼对棒束通道内流场与温度场的影响。

2) 基于搅混翼的阻流、导流、搅混及换热强化等机理所构建的虚拟体积力动量源模型可表示成当地速度场的函数，无需引入其他变量，且数值求解为自动耦合过程，当计算收敛时，合适的动量源大小即被确定。

3) 虚拟体积力动量源模型所计算的横向流场与实验数据及体贴网格模型结果相比，在搅混翼格架近下游区域与二者均吻合较好，远下游区域则出现流速衰减偏大的现象。针对搅混翼棒束通道的压降、Nu、流速等关键热工水力参数，虚拟体积力动量源模型计算结果在数值与分布趋势上均与贴体网格模型结果基本吻合，局部最大偏差在5%以内。此外，虚拟体积力动量源模型对棒束通道来流速度不敏感，能适应不同流速工况。验证结果证明了虚拟体积力动量源法的合理性与可行性。

本文研究主要针对搅混翼，后续将研究其他复杂结构(如格架带条、刚凸与弹簧等)的类似建模方法，有望将燃料组件简化为裸棒束通道而不丢失搅混翼格架的主要热工水力特性，进一步提高精细化建模与仿真效率。