铅基快堆SGTR事故下热工水力模拟及气腔扩散行为研究

于启帆，赵亚峰，王成龙，张大林，田文喜，秋穗正，苏光辉

(西安交通大学 能源与动力工程学院，陕西 西安 710049)

美国能源部(DOE)于2002年与其他12个国家共同成立了第4代国际论坛(GIF)，旨在维持核能的发展。GIF共推出6种四代反应堆堆型，其中铅基快堆(LFR)被认为是最具运行可靠性和性能稳定性的反应堆[1]。考虑到蒸汽发生器传热管两侧存在较大的压差，以及反应堆一回路长期运行造成的铅基合金磨蚀传热管壁面，铅基快堆存在发生蒸汽发生器传热管破裂(SGTR)事故的潜在风险。

针对铅基快堆SGTR的实验研究主要为整体性事故模拟实验研究。意大利ENEA-Brasimone实验室先后通过LIFUS5和LIFUS5/MOD2实验装置，研究了铅铋-水相互作用SGTR过程中的铅池场压力攀升曲线，以及气泡经过热电偶时带来的温度变化[2-4]，此外还搭建了CIRCE系列实验装置，探究了第4代铅铋反应堆MYRRHA的等比例一次侧换热器(PHX)发生SGTR事故时的实验现象[5-6]。德国卡尔斯鲁厄理工学院通过KALLSTARR实验项目，实现欧洲工业嬗变设施(EFIT)1∶4缩比蒸汽发生器的SGTR事故和冷却剂-冷却剂相互作用(CCI)蒸汽爆炸计算[7]。俄罗斯诺夫哥罗德国立技术大学(NNTU)连同俄罗斯物理动力工程研究院(IPPE)搭建了FT-216实验平台，探究了含液滴与不含液滴气泡在铅池内上升到铅池自由液面过程中对铅池的降温影响，进而确认了气泡内水的导热系数[8]。上述实验台架均为SGTR事故全阶段的整体性模拟，不能很好地对事故过程产生的四类热工水力现象进行剥离分析。

本文讨论LIFUS5/MOD2实验台架的SGTR实验工况，建立基于计算流体力学(CFD)软件的两相流模型，通过欧拉-欧拉方法对实验产生的压力波和压力攀升过程进行准确模拟，并对若干不同SGTR事故工况下的压力攀升规律、铅铋合金冲刷扩散规律和水蒸气分布规律进行探讨。

1 意大利LIFUS5/MOD2实验装置

LIFUS5/MOD2实验装置旨在研究欧洲第4代铅基快堆的热工水力问题，实验通过注射管线将高压射流水注射进憋压的反应容器铅池内，观测容器内的温度场和压力波动。反应容器内部的测点布置如图1所示，其中，容器内部共布置4个压力传感器，用于测量高压射流水对反应容器造成的瞬时升压，包括3个高位压力传感器布置于反应容器的东、西、南3个方向的容器内壁面上，1个低位压力传感器布置于反应容器的北方向的容器内壁面上。此外，容器内部另有布置36个热电偶，用于判断低温的水蒸气在铅池内部上升过程的路径，以及蒸汽腔的升温情况。热电偶布置架共计4层，为十字水平形状，中心轴与反应容器同轴，其中自下而上的第1层中心不布置热电偶，其余每层热电偶布置架的分支上各布置4个热电偶。反应容器及实验工况主要设计参数列于表1。

图1 反应容器内部测点布置Fig.1 Lay-out of measuring point in reaction vessel

表1 反应容器及实验工况主要设计参数Table 1 Design parameter of reaction vessel and test case

2 数值模拟模型设置

研究对象为LIFUS5/MOD2实验台架反应容器，通过ANSYS ICEM CFD商业软件完成建模。考虑到LIFUS5/MOD2实验台架反应容器为器身圆柱形、器底半球形的压力容器，满足几何对称性，因此使用一个沿径向贯穿压力容器中轴线的二维切片开展数值模拟研究，并进行热工水力学分析。

2.1 计算域

反应容器模型包含圆柱形器身、缺心的半球形器底和插入容器底部中心的细长注射段，ANSYS ICEM CFD可对上述二维模型生成统一的以压力容器中轴线为轴的轴对称结构化网格，如图2所示。采用了3种不同数量的网格方案进行网格敏感性分析，主要对比参数为：1.6×106Pa的背压水注射开始后前0.1 s的反应容器内最大压力峰值、第0.1 s的升压压力，以及二者对比初始压力的变化百分比。网格敏感性分析列于表2，最终选取157万网格作为本次计算模型的网格规模。

图2 反应容器结构化网格Fig.2 Structured mesh on reaction vessel

2.2 铅铋合金热物性

铅铋合金的热物性对于SGTR事故的热工水力分析具有重要意义。铅铋合金作为铅基快堆的设计冷却剂，由44.5%的Pb和55.5%的Bi组成，熔点为397 K左右，沸点为1 943 K左右，国际原子能机构(IAEA)对铅铋合金的热物性进行了文献综述[9]。根据现有的实验数据和拟合关系式，得到铅铋合金液态金属热物性。

铅铋合金密度ρL为：

表2 网格敏感性分析Table 2 Mesh sensitivity analysis

ρL=11 096-1.323 6T

(1)

式中，T为液态铅铋合金温度。

铅铋合金导热系数λL为：

λL=3.61+1.517×10-2T-1.741×10-6T2

(2)

铅铋合金热容hL为：

hL=159-2.72×10-2T+7.12×10-6T2

(3)

铅铋合金动力黏度系数ηL为：

ηL=(4.56-7.03×10-3T+

3.61×10-6T2)×10-3

(4)

2.3 条件设置

本文研究采用ANSYS Fluent商业计算流体软件实现SGTR事故过程的模拟，在确认选取的数值计算网格后，即可开展SGTR热工水力分析的条件设置。其中，注射口为入口边界，注水温度固定为473.15 K，注水压力随注射背压变化。考虑到SGTR事故过程为憋压过程，此次数值模拟不研究事故泄压后的现象，因此反应容器壁面和注射段壁面均为绝热无滑移边界条件，事故过程中的热损失忽略不计。

2.4 多相模型设置

研究采用欧拉-欧拉两相流模型，在该模型中，两种或两种以上的非混相流体共享1组动量方程，并在计算域和计算单元中跟踪每种流体的体积分数，因此该模型的一个典型应用是在稳态或瞬态条件下跟踪气液界面，也即适用于模拟SGTR蒸汽腔的产生和浮力运动。研究忽略了对闪蒸现象的直接模拟，而是通过计算每个时间步长下的高压水注射量，并将射流水的闪蒸潜热通过用户自定义函数(UDF)作为吸热源项引入到铅池内，同时向铅池内注射与高压水同等质量和背压的水蒸气，实现闪蒸过程带来的吸热模拟和射流模拟，其中水蒸气采用R-K实际气体状态方程考虑其可压缩性。

1) 表面张力

表面张力是作用于铅铋合金表面并使其表面积缩小的力，其值与温度相关，水蒸气分子间的距离大于铅铋合金液体内部的分子间距离，因此水蒸气分子间的相互作用表现为引力。本文研究采用连续表面张力模型(CSF)求解表面张力[10]，该模型认为表面张力是存在于相界面上的、垂直于相界面的连续三维作用力，且与表面张力系数σ有关。针对铅铋合金与水蒸气的任一界面处的表面张力：

(5)

ρ=αGρG+(1-αG)ρL

(6)

对于铅铋合金表面张力系数σ的计算，Keene[11]认为液态金属的表面张力随温度的升高而降低，且在一定温度范围内呈线性关系：

(7)

式中：σm为液态金属熔点下的表面张力；Tm为液态金属的熔点。

根据Alchagirov等[12]列举的现有熔融铅铋合金表面张力试验研究的汇总，可拟合得到液态铅铋合金在1 400 K以下的表面张力σL为：

σL=0.437-6.6×10-5T

(8)

2) 蒸发模型

在流体体积模型的两相流计算中，Lee模型可用来进行蒸发-冷凝的相间传质传热计算，本文将Lee模型通过用户自定义函数引入。Lee模型的水相蒸发冷凝由水蒸气的输运方程控制：

(9)

(10)

式中：αW为液相水体积分数；ρW为液相水密度；coeff为松弛因子，用于调节液相水和水蒸气的质量交换速率，同时可提高计算的收敛性。coeff与气泡直径dG呈反比关系，但二者的值均由于无法确定，因此需通过经验微调coeff以满足实验结果和模拟要求，一般可调范围在0.1～103。

3 计算结果分析

3.1 基准工况分析对比

在实际操作中，实验初始条件与设计工况存在一定偏差，设计工况与初始条件的对比列于表3。

数值计算以实验的真实初始条件为基准工况进行分析研究，并与ENEA-Brasimone实验室的SIMMER-Ⅲ二维多相安全分析程序的计算结果进行了对比。基准工况的压力攀升曲线如图3a所示，根据实验数据得到，在注射开始的前1.4 s内，由于压差导致的高压水射流和蒸发，使得反应容器以1.09 MPa/s的速度近似线性地快速升压，直至反应容器与水容器达到压力平衡。随着实验系统的压力平衡，高压水注射阶段结束，此时反应容器压力缓慢提升至实验峰值压力1.942 MPa，并缓慢回落至与水容器相同的压力1.83 MPa直至注水阀关闭，上述压力提升和回落的时长共计约0.4 s。可以看到，CFD程序模拟的线性快速升压曲线与实验曲线的斜率相似，且实现了实验系统压力平衡后的缓慢回落现象，最大压力的计算值为1.934 MPa，相对误差为0.37%。图3b对注射开始后0.07 s的压力变化进行了放大，根据实验用高频压力传感器可知，在破口注射瞬间反应容器内会出现瞬时压力波，CFD程序通过铅铋合金的可压缩模型实现了压力波的传递模拟，压力波峰值达到了5×104Pa，与实验结果符合较好。

表3 设计工况与初始条件对比Table 3 Comparison of design condition and initial experimental condition

1.6 MPa工况压力曲线：a——前2 s；b——前0.07 s图3 基准工况的压力变化Fig.3 Pressure change of benchmark condition

3.2 射流水背压的热工水力影响

除基准工况外，数值计算拓展了3种不同的射流水背压工况，具体工况系数列于表4。图4示出水的背压为3×105、6×105、1.4×106Pa 3组工况下压力测点监测到的压力攀升曲线。由图4可见：工况1.2的水背压最低，因此该工况下的反应容器可在注射开始1.2 s后即达到压力平衡，且在压力平衡前后1 s内的升压差最小；工况1.3和工况1.4需要1.6 s甚至更长时间达到压力平衡，这表明SGTR事故注射破口较小的情况下，注水背压的提升会提高反应容器升压的剧烈程度和注射时长。

工况1.2下蒸汽腔压力云图如图5所示，在工况1.2中注射1 s后，压力传感器位置测得压力达到2.8×105Pa，此时蒸汽腔尚未突破铅池液面进入到上方氩气空间中，蒸汽腔内压力为2.3×105Pa左右，压力低于附近铅池压力。该现象是铅铋堆SGTR事故中独有的现象，由于铅铋合金的高密度、大表面张力的物理特性，蒸汽腔的成长和上浮会比在水等其他工质中更慢，同时由于反应容器在SGTR实验过程中处于密闭状态，压力会不断攀升，因此蒸汽腔周围的高压力铅铋会对蒸汽腔进行包裹，阻止蒸汽腔的进一步不稳定扩散。

表4 不同射流水背压工况Table 4 Different condition of water pressure

a——工况1.2；b——工况1.3；c——工况1.4图4 不同工况下测点位置压力的变化Fig.4 Pressure change of different conditions at measuring point

a——体积份额；b——压力图5 蒸汽腔压力云图Fig.5 Pressure contour of steam cavity

3.3 铅池液面波动与气腔扩散运动

通过计算SGTR事故过程中的铅池液面波动，对液面卷吸和被卷吸氩气气体的运动轨迹进行了分析，同时比较了不同射流水背压带来的铅池液面攀升现象，并对各测点空泡份额变化和蒸汽腔的扩散速率差异进行了对比。

1) 铅池液面波动与气体卷吸现象

随着蒸汽腔在铅铋池内的迅速扩大和上浮，铅铋池的中心液面会优先被蒸汽腔挤压和上升，形成中间高、两边低的液面形状，并发生液面卷吸现象。图6示出工况1.2的铅铋液面卷吸气体的过程。液面卷吸现象发生在铅铋池与反应容器的接触位置，即铅铋池液面边缘，这是由于注射口在反应容器中间，导致中间液面的铅铋合金先被冲起，对上方氩气空间造成挤压和对容器顶部造成冲击，并在回落时卷吸边缘液面上方的氩气。氩气腔被卷吸后会跟随铅铋合金的速度惯性而向下俯冲进入铅池，并随后出现两种运动现象：(1) 若被卷吸的氩气靠近蒸汽腔，则会与附近的蒸汽腔合并，并共同上浮至反应容器顶部的氩气空间中；(2) 若被卷吸的氩气距蒸汽腔较远，则会以小气腔或小气泡的状态在铅池内发生夹带运动，并有可能发生气泡尺寸的二次分割。

图6 铅铋液面卷吸气体过程Fig.6 Process of argon entrained by LBE

在工况1.2中被卷吸氩气腔的最小尺寸为10 mm，在更高背压的工况1.4下，被卷吸氩气腔的最小尺寸可达到3 mm。根据Yu等[13]的结论，直径为0.3～0.6 mm和1 mm的汽泡在ELSY铅铋反应堆中具有最大的被夹带至堆芯的概率，这一汽泡尺寸和夹带概率与一回路铅铋合金流动速度和反应堆几何结构相关。

2) 蒸汽腔扩散过程

在LIFUS5/MOD2实验蒸汽腔的扩散和上浮过程中，通过热电偶监测到的温度变化，可近似判断出该测点有水蒸气经过，进而近似反应出蒸汽腔的经过路径。而在本文的数值模拟中，可直接提取到计算域内任意位置的各相体积份额变化，无需通过温度的波动来判断。为保证与实验的一致性，数值模拟选取与实验布置的热电偶相同的位置作为监测体积份额变化的测点。图7示出工况1.3的不同热电偶测点位置，在注射开始后第0.2、0.5、1.0和1.5 s时刻测到的蒸汽空泡份额，其中L1～L4分别表示从低到高的4层热电偶，C、R1、R2、R3、R4分别表示该层热电偶在中心处以及距离中心依次远离的测点。由图7可见：(1) 随着注射时间的推进，从下至上的各体积份额监测点依次监测到水蒸气的空泡份额；(2) 第1层热电偶远离中心位置，在注射开始前1.0 s并未检测到水蒸气空泡份额，而第2、3层的各热电偶在同一时刻测得水蒸气空泡份额的结果更相近，其中第2层热电偶在t=0.5 s时测得的空泡份额在0.38～0.7之间，第3层热电偶在t=1.0 s时测得的空泡份额在0.24～0.55之间，这表明中间高度的蒸汽腔开始沿径向扩散，且扩散均匀；(3) 第4层的C、R1、R2 3个位置的热电偶，在t=1.0 s并未测得水蒸气空泡份额，而R3、R4测到了空泡份额，这表明蒸汽腔上浮到一定高度后，蒸汽腔会沿径向向远离中心的位置扩散，中心处空泡份额降低。

通过以上3点可知，蒸汽腔在注射初始阶段(前1.0 s)呈喇叭形扩散，即轴向高度越高，蒸汽腔在径向的扩散越剧烈，同时随着蒸汽腔的上浮，径向中心处的蒸汽腔空泡份额开始降低。上述扩散规律会随着时间的推移(t=1.5 s)和铅铋液面的卷吸而消失，此后蒸汽腔在铅铋池内的分布趋于均匀。

图7 各层热电偶测点的蒸汽腔空泡份额的变化Fig.7 Volume fraction change of steam cavity at measuring point of thermocouples in each layer

4 结论

本文采用CFD模型对意大利LIFUS5/MOD2实验装置的高压水射流铅铋热工水力SGTR实验进行了研究，采用了铅铋和水的两相模型，并考虑了水的蒸发传质传热计算。将CFD计算结果与实际工况进行了对比研究，主要研究结果如下。

1) 通过基准研究，对于反应容器的压力攀升过程，CFD模型可更好地模拟到注射开始后0.05 s内的5×104Pa的压力波传递现象，且最大压力计算相对误差为0.37%。对于蒸汽腔的扩散与上浮过程，CFD模型存在提前模拟到蒸汽腔上浮的现象，但可更好地预测到测点位置的空泡份额峰值与空泡份额波动现象。

2) 在不同水容器背压引发的SGTR事故中，射流初始时刻的压力波峰值会随水侧背压的升高而增大，水侧背压越高，注射期间的反应容器升压越剧烈。在工况1.2的蒸汽腔上浮过程中，发现蒸汽腔压力(2.3×105Pa)低于附近铅池压力(2.8×105Pa)，这表明压力不断攀升的高压力铅铋会对蒸汽腔进行包裹，阻止蒸汽腔的不稳定扩散，该现象可降低蒸汽腔在铅铋池内出现蒸汽爆炸的可能性。

3) 蒸汽腔在注射初始阶段呈喇叭形扩散，即蒸汽腔沿径向的扩散随轴向高度的增高而剧烈，同时随着蒸汽腔的上浮，径向中心处的蒸汽腔空泡份额会开始向边缘转移，即蒸汽腔上部边缘的空泡份额高于中间。该扩散现象可持续1 s，并会随时间的推移而消失，此后蒸汽腔在铅铋池内的分布趋于均匀。

本文研究未考虑压力波的传递对结构件造成的力学冲击和形变分析，后续研究将进行进一步的深入探索。