燃料组件堵流工况下铅铋-氩气两相流的传热压降特性分析

2024-03-01刘书勇李桃生梅华平赵吉运

核安全 2024年1期

夏凡，刘书勇，李桃生，梅华平，汪振，赵吉运

（1. 中国科学院合肥物质科学研究院核能安全技术研究所，合肥 230031；2. 中国科学技术大学，合肥 230026；3. 香港城市大学机械工程系，香港）

第四代铅冷快堆具备安全性好及易小型化的优势，在分布式供电、海岛平台等领域具有广阔的应用前景。在安全性方面，无论是回路式还是池式设计，堵流都是重要的事故工况，必须给予考虑［1-3］。堵流的起因主要是由于氧控不良引起氧化铅在组件内的沉积［4］，或是由于腐蚀脱落的结构材料碎片随冷却剂带入，造成组件内的堵流。燃料组件内堵流发生后会进而引发过热和局部流量下降，影响堆的安全运行［5-8］。

对于液态铅铋冷却的堵流实验研究，Pacio J 等开展了块状堵块实验，通过设置低导热系数的固体堵块，研究了不同堵流工况下的流动传热［9］。该堵流实验指出，实际堆运行中堵块多为多孔介质，由于多孔堵块的加工制造存在困难，因此该实验采用的堵块并没有设置孔隙率。2019 年，意大利ENEA 基于铅铋冷却回路（NACIE-UP）开展了堵流实验研究，结果表明，堵流后方形成局部温度峰值，堵流子通道活性区末端形成全局温度峰值［10］。对于堵流模拟研究，2014 年，Di Piazza 等建立CFD 模型对无绕丝工况下燃料组件堵流现象进行了数值模拟，研究结果表明：对于多子通道堵流，在堵块下游的回流区域出现温度峰值；对于少子通道堵流，温度峰值出现在活性区域末端［11］。2019 年，上海交通大学柴翔等基于KIT 实验装置，采用多子通道堵流方式进行了模拟，结果表明中心通道堵流包壳温度高于角通道温度［12］。2020 年，吕科锋等对带绕丝19 棒束铅铋多孔介质堵流进行了模拟，得出包壳存在圆周温度变化［13］。上述研究均限于铅铋冷却单相流，未涉及两相流的流场计算。

对于铅基堆而言，通常采用氩气作为保护气体或在加强自然循环中提供驱动压力［14-16］，因此存在铅/铅铋-氩气两相流工况［17-19］。本文对于液态金属-氩气两相流的数值模拟主要关注如下两个方面的研究。一方面，对于加强自然循环工况［20-23］，关注的重点在于铅铋空泡份额与氩气注入速率之间的关系［21，23］。另一方面，在泵驱动的强迫循环工况下，回路运行前需先用氩气将空气排空，再注入氩气将熔化的液态铅铋打压进入回路。组件复杂结构内的部分区域难以被铅铋充满，残留的氩气会滞留在回路内并形成局部“死区”，可能因导热不良引发局部过热，如有堵流发生，可能会加剧局部温升。

因此，本文通过数值模拟给出铅铋-氩气两相下的气泡行为，进而分析堵区局部过热和压降特征。

1 研究方法

1.1 单相流计算方法

通常，燃料组件内液态铅铋的单相流动可以视为不可压缩流体，满足连续性、动量守恒和能量守恒方程。文献［33］针对采用气泡提升泵加强自然循环工况，对比了ε型及ω型湍流模型的相对误差，得出standardk-ε计算误差较小且节省计算资源，因此本文采用standardk-ε湍流模型，详细公式可以参见ANSYS FLUENT theory guide［24］。对于多孔介质堵流，堵塞区域要增加考虑动量源项，见式（1）～式（3）。其中，Sporous是多孔区域动量源项。

式中，下标hetr.可以是i，j，k三个方向的分量，表示速度在空间x，y，z上呈各向异性变化。多孔区域内的动量损失由黏性损失［即式（1）第一项］和惯性损失［即式（1）第二项］构成［25］。式（1）中变量β为不透水性，定义见式（2）。

式中，dpar是颗粒直径，mm。是孔隙率。该方法首次由Ergun 于1952 年提出［26］，并且可以应用于铅铋多孔介质模拟仿真［17］。式（1）中惯性损失常数定义为摩擦系数和水利直径的比值，如公式（3）所示。

针对液态铅铋低普朗特数特性，即液态铅铋与水相比具有更大的热边界层［27］。在ANSYS FLUENT 软件中，默认的对流换热和湍流普朗特数是针对液态水进行求解的，因此在求解铅铋介质时会带来误差。Lyu Kefeng等针对P/D为1.14、带绕丝的铅铋冷却燃料组件堵流工况，采用如公式（4）［28］所示的湍流普朗特数（Prt）关系式，获得了较为合理的分析结果［13］。

1.2 两相流计算方法

两相流的处理方法分为两类：一类是均相流模型，适用于气相含量比较低，并且两相相对速度不大的情况；另一类是分相流模型，即对气液两相分别进行处理，对每一相计算其平均物理参量。ANSYS FLUENT 提供的多相流处理方法有VOF（Volume of Fraction）模型、混合物（Mixture）模型和欧拉（Eulerian）模型，前两种模型属于均相流模型，而Eulerian 模型属于分相流模型。鉴于泵驱动的强迫循环和加强自然循环工况下，氩气含量很少，铅铋-氩气相对速度不大，因此可以采用VOF 模型，采用CSF（Continuum Surface Force）方法模拟表面张力［29］。

2 模型建立及验证

2.1 铅铋冷却燃料组件单相流模型

本文以德国KIT 铅铋冷却燃料组件实验流动传热数据为基准［9，32］，进行数值计算结果的比对验证。实验装置侧视图及堵块设置如图1 所示，燃料组件主要参数见表1。堵流实验采用电加热的均匀面热源，热功率为394±4 kW。铅铋入口温度为200±0.2℃，质量流量为18.7±0.2 kg·s-1［32］。

图1 燃料组件实验装置侧视图（黄色矩形区域为堵块位置）［11］Fig.1 Side view of the testing platform （blockages are highlighted in yellow）［11］

表1 燃料组件主要参数［32］Table 1 Main parameters of the fuel assembly［32］

本课题组已完成了正常运行工况和全堵工况的流动传热及网格无关性验证，结果显示，温度的模拟结果与实验相比，最大相对误差在±4.9%，摩擦系数与实验推荐经验关系式的相对误差在±9.2%［33］。相对误差的计算方法为：（数值模型计算值-实验参考值）/实验参考值。

2.2 铅铋-氩气两相流模型

基于表1 的燃料组件结构，本文采用FLUENT 软件中自带的FLUENT Meshing 模块，针对表1 的中心6 子通道堵流（简称C6工况）和边1 子通道堵流（简称E1工况）两种工况，分别进行几何结构建模及网格划分，并在堵块上游设置一个圆形气泡的起始位置，如图2所示。

图2 堵块及气泡起始位置局部网格划分Fig.2 Mesh display for the blockage regions and gas bubbles

加热棒的壁面边界条件采用与图1 基准实验一致的均匀热流密度，对于小堵块E1工况，为0.93 MW·m-2；对于大堵块C6工况，考虑降功率保护，设置为E1工况的1/4，即0.23 MW·m-2。铅铋入口条件采用不同雷诺数下的质量流量进口，温度为473.15 K；出口处设置为压力出口，压强取0 Pa。氩气气泡的起始位置设置在堵块区域前方10 mm 处，直径为1.0～2.0 mm（见表3 第2 列），并对气泡可能流经的区域进行网格加密。本文采用FLUENT 软件对流场进行初始化，在圆形气泡区域通过patch 方式使其充满氩气。在瞬态计算时，设置连续性方程收敛条件取10-4量级，动量及能量方程收敛条件取10-5量级。考虑到流速和最小网格尺寸，时间步长取1×10-3s。

为验证两相流模型，本文同步通过速度拟合法［34］和气泡稳定形态图示法［16］进行验证。其中，前者属于经验关系式，后者属于经验图示，均来源于对大量实验数据和结果的总结归纳。而本文数值模型采用VOF 方法属于两相求解方法，其与连续性方程、动量方程和能量方程一并用于两相速度、能量的计算求解。采用VOF 方法得到的数值模拟结果若与经验图示及经验关系式的求解结果相一致，则表明VOF 方法用于铅铋-氩气工况计算是适合的。

对于速度拟合法，研究指出，通过计算无量纲数Fr，Eo，Mo和Nf数（Nf数为Eo、Mo的组合）可以得到两相间作用力的影响情况。当Eo＞70 时，表面张力可以不考虑；当Nf＞550时，黏性作用影响很小；当Fr＜0.05 时，惯性作用可以不考虑［34］。根据ANSYS 理论手册［24］，对于Re＞＞1 的工况，We数（惯性力和张力之比）不可忽略，当We＞＞1 时，张力影响可以忽略。

这些无量纲数的表达式（见表2 第2 行）与本文工况点1～6（见图4）的计算结果各参数列于表2 第3 行。由计算结果可知，惯性力和张力的影响是主要影响因素，惯性力占主导。当惯性力的影响占主导时，根据文献［34］可知，垂直管内流体速度应当满足u=Const.×Ul+Ug，c，其中Ul为液相流速，m/s；Ug，c为最终稳定时气泡头部中心线上的速度，m/s。对于湍流，Const.的经验值取1.2。由于VOF 模型只能获得流体速度u，因此，将流体速度u与铅铋入口速度Ul进行拟合（见图3），得到Const.约为1.05，该值稍小于经验值［34］，相关系数R2为0.9997，这与其他研究者的结论一致，因此可以认为模拟的结果是可以接受的。

图3 惯性力主导下流体速度和液相速度拟合结果Fig. 3 The fitting curve for inertial force dominated velocities between two-phase flow and the liquid LBE

图4 本模型工况点1～6 气泡稳定上升时的形态Fig.`4 The numerical results of rising bubbles for E1 blockage conditions

表2 无量纲数的表达式及本模型中计算结果Table 2 Expressions for dimensionless numbers and calculation results

对于气泡稳定形态图示法［30，31］，本文选取边1 子通道堵流工况下的6 个工况点作为研究对象，其中1～2 研究气泡直径的影响，3～6 研究雷诺数的影响，工况点各参数取值见表3。网格加密区域的两相流流速和气相雷诺数由FLUENT 求解器输出，分别见表3 第4 列及第5列。Mo数、Eo数、Re数定义式分别见式（5）～式（7）［16］，通过计算获取这三个无量纲数列于表3后3 列之中。由FLUENT 软件计算后，获得气泡随无量纲数变化稳定时的形态如图4 所示。考虑到计算结果的可靠性，选取Grace 经验图示（见图5）作为基准形态，将表3 参数对照图5 坐标，画出工况点1～6 各点位置，见图5 黄线所示。根据图5 形态划分，工况点1 落在球状（Spherical）区域，其余5 个工况点均落在抖动状（Wobbling）区域。对比图4 中本模型中输出的稳定形态，以及图5 中气泡稳定时的理论形态，可知本论文采用VOF 方法的模拟结果与经验结果相吻合。

图5 Grace 经验图示下的理论形态Fig.5 The stable shapes for gas bubbles under Grace correlation

表3 氩气气泡在液态铅铋中稳定上升时各参数取值（E1 工况下，孔隙率0.2）Table 3 Values of parameters for steady rising bubbles for E1 blockage

3 堵流工况下的两相流特性分析

本节考虑铅铋夹带氩气工况下，基于两相流模型验证结果，对C6工况和E1工况的局部过热和压降分别进行计算分析，并对入口雷诺数、孔隙率及气泡直径进行了参数敏感性分析。堵块边界条件设置分为两种：一种是以外边界为壁面来模拟实心堵块情况；另一种是外边界为内部面并设置孔隙率来模拟多孔介质堵块情况。

3.1 C6 工况

3.1.1 入口雷诺数的影响

（1）堵块外边界为壁面条件

本文通过设置气相等体积分率的监测面，观察含气区域因导热不良造成的温升。选取组件轴向位置348～552 mm 的体平均温度为参考温度，堵块在492～546.6 mm 区域，相对于参考节段位置为70.6%～97.4%，气泡直径保持2 mm不变。图6 给出了不同雷诺数下气相等体积分率（VOFAr）在［0.6～0.9］范围内的过热情况。由图6 可知，随流动时间的增加，VOFAr较高的等值面不一定存在连续性。其中图6（a）至图6（b）表明雷诺数较低时，VOFAr为0.9 的等值面大概率能够在0.08 s 内监测到；而雷诺数较高时，且VOFAr值较大时，气泡存在时间较短，仅约为0.03 s，如图6（c）所示。且VOFAr值越高，引起的局部过热越明显。壁面边界与内部面边界相比，在瞬态条件下气相体积分率可以保持在较高的水平（CAr＞0.5），从而导致引起的过热更严重。

图6 不同雷诺数下气相等值面温度变化（C6 工况，外边界为壁面）Fig.6 Temperature of iso-surfaces for gas phase at C6 blockage conditions（wall boundary conditions）

图7 给出了当流动时间为0.01 s，Re为15000 时入口条件下的局部VOFAr分布及温度云图。由图6 的计算结果可知，尽管VOFAr监测面温度均值约在550 K，尚在可接受的范围，然而图7 局部最高温度接近1000～1500 K，这可能对局部结构完整性造成威胁。

图7 Re =15000 时堵区局部过热，t=0.01s （C6 工况）Fig.7 Local overheating for diff. Re of C6 blockage conditions

图8 给出了中低雷诺数下气泡被封闭的现象，由图可知雷诺数较低时，气泡可能会卡在堵块附近形成“死区”，并伴随明显的过热。图9 显示了高雷诺数下气泡逸出现象，由于在高雷诺数下高速流动的液态铅铋可能将气泡带出堵流区域，引起气泡逃逸，从而使高含气量的维持时间很短。即便在0.1 s 后再次监测到高等值面，逸出气体引起的局部过热也微乎其微。

图8 中低雷诺数下气泡被封闭现象（C6 工况）Fig.8 The sealed bubble for low and immediate Re of C6 blockage conditions

图9 高雷诺数下气泡逸出现象（C6 工况）Fig.9 The escape bubbles for high Re of C6 blockage conditions

（2）堵块外边界为内部面条件

在该类边界条件下，气泡可以贯穿多孔堵流区域，最终逃逸出该区域。为了研究不同入口Re影响，本文取孔隙率为0.8 和气泡直径为2 mm 且保持不变，观察不同时刻气相等值面温度变化。

图10 给出了不同雷诺数下VOFAr等值面温度计算结果，由图示可知，高VOFAr值面存留的时间非常短，且低VOFAr值时温度随流动时间近似线性升高。例如，Re取5000 的条件下，VOFAr取0.8 的等值面存留时间仅为50 ms，如图10（a）所示。由图10（b）可以看出，随着雷诺数增加，VOFAr等值面存在时间缩短，且最终引起的温升有所降低。然而，堵块在内部面边界条件下引起的过热并不显著。例如，流动时间为0.15 s 时，以堵块所在节段的流体的体平均温度为参考，在Re为5000 时VOFAr取0.6 的等值面过热为5.4 K，而壁面条件的相同雷诺数下，其过热温度已达22.6 K。

图10 不同雷诺数下VOFAr 等值面温度计算结果（C6 工况，外边界为内部面）Fig.10 Temperature of iso-surfaces for gas phase at C6 blockage conditions （internal boundary conditions）

为监测氩气及堵流区域的静压瞬态变化，研究者在FLUENT 软件中创建了一条监测线贯穿多孔堵流区域。该直线的x，y坐标与气泡起始位置的坐标保持一致，得到不同雷诺数下沿轴向方向静压变化曲线。图11 显示随着堵块位置轴向坐标的增加，高VOFAr区域可能引起局部静压的小幅上升（微正压现象）。特别是在气泡进入多孔区域瞬间，引起的相对静压增量约为0.3%。由图11（a）～图11（b）的变化可知，雷诺数越低，微正压现象越明显。这一现象从两相流稳定条件出发可以解释为：气泡在液体中平衡存在除了需要具有一定的过热度外，气相和液相压差还需满足pg-pl=2σ/r*条件，其中pg为气泡内的压力，Pa；pl为液相压力，Pa；σ为表面张力，N/m；r*是界面曲率，1/m［35］。因此，对于不透明介质，可以利用这一现象来监测气泡的位置。

图11 不同雷诺数下两相区域压降计算结果，Re 5000～15000（C6 工况）Fig.11 The pressure drops of diff. inlet Re for C6 blockage conditions

3.1.2 孔隙率的影响

为研究孔隙率对局部过热和压降特征的影响，本文铅铋入口雷诺数取临界Re，即过渡流向湍流转变的RebT=15700，且气泡直径取2 mm 保持不变，计算得到不同时刻孔隙率［0.2，0.8］的VOFAr值监测面温度曲线，如图12 所示。

图12 不同孔隙率下温度计算结果（C6 工况，Re=15700，气泡直径为2.0 mm）Fig.12 Temperature of iso-surfaces for gas phase at C6 blockage conditions

由图12（a）可知，孔隙率取0.2 时，等值面温度比该节段流体的体平均温度过热显著，温度曲线和全堵工况趋势类似。这是由于气泡被限制在堵块局部形成 “死区”［见图13（a）］，与全堵工况时形成的“死区”相比，低孔隙率下气泡并未完全被限制，而是分裂成一个较大的主气泡和一个较小的子气泡。主气泡未被限制而逃逸，而子气泡被限制在堵块区域，如图13（b）所示。

图13 孔隙率0.2 时气泡行为，t =0.1 s（C6 工况）Fig.13 The behavior of gas bubbles when porosity is 0.2 at C6 blockage conditions，t =0.1 s

对于孔隙率为0.4［见图12（b）］的情况，在0.03 s 后仅低VOFAr值能够被监测到，这种工况下，气相也会引起比较明显的过热。气相随流动时间增加分裂成多个子泡，缓慢逃逸出多孔堵流区域。当孔隙率增至0.6 时，除了可以观察到气泡逃逸［图14（a）］之外，气相还会随流动时间增加缓慢耗散在多孔堵流内部，如图14（b）所示。当孔隙率增至0.8 时，气泡能够完整地随流体从堵块出口逸出多孔介质区域，并不会引起明显的过热现象。

图14 孔隙率0.6 时气泡逃逸和耗散现象（C6 工况）Fig.14 The behavior of gas phase when porosity is 0.6 at C6 blockage conditions

3.1.3 气泡直径的影响

为探究气泡直径的影响，研究者分别设置直径为2.0 mm、1.5 mm 和1.0 mm 的圆形气泡，其余参数取Re=RebT=15700，由于小孔隙率可能引起较为严重的过热，故选取孔隙率为0.2 并对所有工况保持一致。

图15 为不同气泡直径下VOFAr值监测面温度计算结果。由图15（a）可知，仅在直径为2.0 mm 时，高VOFAr值监测面的存留时间较长，在0.02 s 内引起的过热较为明显。图15（b）显示直径为1.5 mm 的气泡在0.015 s 左右存在约为2 K 的短期过热。这是由于仅直径为2.0 mm 分裂出子泡并伴随有“死区”形成（见图13），其他气泡直径下均未分裂出子泡。因此，在铅铋-氩气的夹带工况中，应当关注气泡尺寸（直径≥2.0 mm）较大的情况；在氩气加强自然循环工况中应当合理控制气泡的注入尺寸。

图15 不同气泡直径下气相温度计算结果（C6 工况，孔隙率0.2，Re=15700）Fig.15 Temperature of gas iso-surfaces for different bubble diameters（ =0.2，Re =15700，C6 blockage conditions）

3.2 E1 工况

3.2.1 入口雷诺数影响

与上述探讨C6工况时控制变量法的思想类似，为研究铅铋入口雷诺数影响，研究者控制气泡直径为2.0 mm，孔隙率为0.8，设置入口雷诺数以5000 为增量递增，变化范围为5000～40000，观察堵流工况下气泡行为、压降及过热特征，结果如下。

（1）堵块外边界为壁面条件

在边1 子通道堵流（E1工况）下，气泡相较中心通道堵流有更大的逃逸趋势。图16（a）（b）分别显示了Re为10000 和15000 时气泡的逃逸现象。模拟结果显示，在Re≥10000 时，气体在多孔区域内的存留时长不超过0.01 s 便缓慢逸出。因此，除了Re为5000 时气相的过热较为明显外，其余雷诺数下气相引起的过热皆不明显。

图16 不同雷诺数下气泡逃逸现象（E1 工况，外边界为壁面）Fig.16 Different behavior of gas bubbles for different inlet Re of E1 blockage

（2）堵块外边界为内部面条件

在E1工况下，堵块位于轴向710.7～765.3 mm段，参考节段为组件加热段轴向696～870 mm时，堵塞位置所在该段8.4%～39.8%的中上游位置。图17 显示了Re为20000 的条件下，堵块区域静压的变化情况。由图示可知，当流动时间由0.01 s 增至0.03 s 时，高气相含量区域存在与中心6 子通道工况类似的微正压现象（图中箭头所指位置为气泡位置）。

图17 Re=20000 时静压变化情况（E1 工况，外边界为内部面）Fig.17 Pressure changes for E1 blockage，Re=20000

温升方面，在0.03 s 内，VOFAr等值面温度近似指数上升，然而其引起的过热并不明显，因此本文不进行展开阐述。

3.2.2 孔隙率影响

边通道堵流下，孔隙率大小可能会对气泡在多孔介质内的存留时长造成影响。在图18中，研究者对不同孔隙率下气泡是否逃逸进行了对比，由模拟结果可知，当孔隙率不太大（＜0.8）时，气泡容易在短期内逃逸，VOFAr等值面温度上升到一个峰值温度然后下降；当孔隙率比较大（ ≥0.8）时，气泡不易从多孔介质内逸出，因此能够被加热较长的时间，VOFAr等值面温度近似指数形式上升，温升也比小孔隙率时更为显著。

图18 不同孔隙率气泡行为，Re=RebT，气泡直径为2.0 mmFig.18 Different behavior of gas bubbles for different porosities，diameter 2.0 mm

3.2.3 气泡直径影响

流动时间为0.01 s 时，气泡直径分别为1.0 mm、1.5 mm 和2.0 mm 的气相参数和两相流动的模拟结果已经在模型验证部分进行了阐述（见图4）。对于不同直径，气泡皆能够存留在多孔区域较长时间，但引起的温升均在1 K 以内，并不显著。随流动时间的增加，不同直径的气泡缓慢耗散在液态铅铋中，VOFAr监测值也逐渐降低。