铅铋堆蒸汽发生器传热管破裂事故下铅铋-水相互作用程序开发及验证

2023-08-01辜峙钘余红星黄代顺严明宇申亚欧张牧昊

原子能科学技术 2023年7期

辜峙钘,余红星,黄代顺,严明宇,申亚欧,张牧昊

(1.中国核动力研究设计院核反应堆系统设计技术重点实验室,四川成都 610213;2.成都理工大学核技术与自动化工程学院,四川成都 610059)

铅铋反应堆对我国核电可持续、安全发展具有重要战略意义[1]。蒸汽发生器传热管破裂(SGTR)事故是铅铋反应堆设计必须重点考虑的关键安全问题之一。蒸汽发生器(SG)是反应堆内承担热传输的核心部件之一,其传热管壁长期处于恶劣环境中,发生破裂频率较高[2]。同时,与压水反应堆相比,铅铋反应堆SGTR事故具有一定特殊性与挑战性,铅铋反应堆的主回路系统通常采用一体化池式结构设计方式[3],其SG、主泵等关键设备直接浸泡在铅铋池内。传热管一旦破裂,高压过冷水将喷射进入一回路内,与高温液态铅铋直接接触,引发铅铋-水相互作用,可能引发一系列安全相关问题[4]。因此从安全角度,SGTR事故是铅铋反应堆设计必须重点考虑的安全问题之一。

铅铋反应堆SGTR事故引发的安全问题,究其根源,在于铅铋-水相互作用行为。铅铋-水相互作用行为是认识铅铋反应堆SGTR事故演化之根本,也是开展铅铋反应堆SGTR事故安全评价的方法学基础。对于铅铋反应堆SGTR事故铅铋-水相互作用行为问题,国内外已开展了一定实验和数值模拟研究。对于实验:一方面,日本原子能研究机构[5-6]、东京工业大学[7-8]、中国科学院[9]等主要开展了铅铋-水相互作用机理性实验研究,主要聚焦于空穴的形成、演化,液态金属碎化,蒸汽爆炸等现象、行为、机理;另一方面,欧洲核能机构、中国科学院等针对铅铋反应堆SGTR事故行为,分别先后搭建了LIFUS5[10-12]、KYLIN-II-S等系列铅铋-水相互作用工程试验平台,具体研究了高压过冷水注入高温、常压液态铅铋之后的宏观演化过程[10-12]。然而,对于铅铋反应堆SGTR数值模拟研究,传统两流体程序并不适用,且针对性程序更是凤毛麟角,国外主要采用液态金属反应堆严重事故分析程序SIMMER-Ⅲ[13],国内主要聚焦于多相流程序NTC-2D[14-15]、MC3D[16],并开展了相关数值模拟研究,但总体而言,相关研究仍处于探究阶段。针对该问题,本文开展铅铋反应堆SGTR事故下池内铅铋-水相互作用理论、关键模型及相关算法研究,开发完成一款铅铋-水相互作用三流体程序,并通过使用国外已公开报道的实验数据进行程序验证。

1 理论、关键模型及算法

1.1 铅铋-水相互作用理论

液态铅铋、液态水、蒸汽混合物(包括铅铋蒸汽、水蒸气、不可凝气体)的质量守恒方程如式(1)～(5)所示,本文采用二维轴对称圆柱坐标系,考虑径向和轴向两个方向。方程右侧考虑了相变传质行为,同时可看出3种蒸汽组分共享同一速度场。

(1)

(2)

(3)

(4)

(5)

对于动量守恒方程,考虑了3个速度场,即液态铅铋、液态水、蒸汽混合物,其径向动量守恒方程分别如式(6)～(8)所示。方程右侧考虑了相间拽曳力、虚拟质量力及相变传质带来的动量交换。轴向动量方程与径向相似,限于篇幅,这里不再赘述。

FLBE,G,R+FVM,LBE,R+ΓC,LBEuG-ΓV,LBEuLBE

(6)

FH2O,G,R+FVM,H2O,R+ΓC,H2OuG-ΓV,H2OuH2O

(7)

ΓC,H2OuG-ΓC,LBEuG+ΓV,H2OuH2O+ΓV,LBEuLBE

(8)

对于能量守恒方程,考虑3种能量组分,即液态铅铋、液态水、蒸汽混合物,它们的能量守恒方程分别如(9)～(11)所示,其中忽略了网格间热传导以及拽曳力、虚拟质量力做功。

(9)

(10)

QΓ,G+QG(h,a,ΔT)

(11)

式中:eLBE、eH2O、eG分别为液态铅铋、液态水、蒸汽混合物的比内能;QΓ,LBE、QΓ,H2O、QΓ,G分别为由于交界面相变对液态铅铋、液态水、蒸汽混合物造成的传热速率;QLBE、QH2O、QG分别为由于交界面温差传热对液态铅铋、液态水、蒸汽混合物造成的传热速率;h为传热系数;a为交界面积;ΔT为温差。

1.2 铅铋-水相互作用关键模型

1.2.1多相流流型模型就本文而言,主要还是考虑池内铅铋-水相互作用行为,并未考虑堆内复杂结构(如蒸汽发生器传热管道、堆芯棒束通道等)环境,同时考虑到铅铋-水接触引发急剧传热、传质,发生泡状流、弥散流、过渡流3种流型演变的可能性较大,因此便考虑了泡状流、弥散流、过渡流这3种典型流型。上述3种流型通过空泡体积份额判断,当αG≤0.3时为纯泡状流,当αG≥0.7时为纯弥散流,当0.3<αG<0.7时为过渡流。对于任意控制体网格,均可看作由泡状流区域、弥散流区域构成,只是上述泡状流、弥散流为两种特殊情况,作为描述流型特征的关键参数,即泡状流区域、弥散流区域体积份额fB、fD按式(12)计算。

fD=1-fB

(12)

如图1所示,为模拟铅铋-水之间急剧传热、传质及水的射流冲击造成的相界面突变问题,考虑了对泡状流区域进行二次划分,考虑第一泡状流和第二泡状流两区域,即液态铅铋、液态水两种液体均可同时作为连续相处理。作为泡状流区域二次划分的关键参数,第一泡状流区域、第二泡状流区域体积份额fB1、fB2按式(13)计算,其中αCP1,e为第一连续相体积份额,对于铅铋-水相互作用而言,若第一连续相为液态铅铋,则αCP1,e=αLBE/(αLBE+αH2O)。

白色连续填充为气体,蓝色与灰色连续填充为液体(铅铋或水),白色圈为气泡,蓝色与灰色圈为液滴图1 多相流流型划分方案Fig.1 Multiphase flow pattern division scheme

fB2=1-fB1

(13)

1.2.2相界面输运模型为模拟在流动、传热、传质等过程中相界面的演化行为,程序求解了相界面输运方程模型,具体包括泡状流区域内气泡和液滴、弥散流区域内液滴三类相界面输运情况,如式(14)～(16)所示。

(14)

(15)

(16)

式中:aG,B为泡状流区域内气泡相界面积;VG为蒸汽混合物速度矢量;αG,B为泡状流区域内气泡体积份额;SG,B,N为泡状流区域内气泡相界面气泡成核源项;aD,B为泡状流区域内液滴相界面积;aD,D为弥散流区域内液滴相界面积;SD,D,Fl为弥散流区域内液滴相界面闪蒸源项。

(17)

(18)

式中:Mb为模化系数,与液体温度、液体饱和温度、液体临界温度、流型参数等相关;pS,D为液滴饱和压力;re,Fl为液滴闪蒸平衡态半径;pD,G为液滴相应的蒸汽分压;σD为液滴表面张力;τD,D,Fl为弥散流区域内液滴闪蒸时间常数;rD,D为弥散流区域内液滴半径;evap,D为液滴汽化比内能;econ,D为液滴相应蒸汽凝结比内能;υvap,D、υcon,D分别为汽化、凝结比体积;λD为液滴热导率;TD为液滴温度;TS,D为液滴饱和温度。

1.2.3相间拽曳力模型相间拽曳力首先根据不同流型条件、不同接触方式进行计算,然后通过加权平均获得液态铅铋、液态水、蒸汽混合物3种流体之间相间拽曳力模型。根据不同的流型条件,接触方式包含连续液体包裹液滴、连续液体包裹蒸汽泡、连续蒸汽包裹液滴、液滴-液滴。据此可分为连续相-弥散相、弥散相-弥散相两种类别,这里以径向拽曳力为例进行说明。

对于连续相-弥散相接触类别:

FCD,R=KCD,R(uC-uD)=

(19)

式中:FCD,R为连续相对弥散相施加的径向拽曳力;KCD,R为其相应的拽曳力系数;uC、uD分别为连续相、弥散相径向速度分量;aCD为其交界面积;μC为连续相动力黏度;rD为弥散相半径;CD为系数,可调;ρC为连续相理论密度。

该模型考虑了黏性力效应和惯性力效应,其中,系数CD参考Ishii等[18]发展的相关模型。

对于弥散相-弥散相接触类别:

FD1D2,R=KD1D2,R(uD1-uD2)=

(20)

弥散项-弥散项接触下,忽略了其黏性效应。轴向拽曳力系数模型与径向相似,本文限于篇幅不再赘述。

1.2.4虚拟质量力模型在两相流或多相流中,当弥散相(如气泡或液滴)相对连续相加速时,将产生虚拟质量效应,动量方程中的虚拟质量力便是为了描述该效应。而且在铅铋堆SGTR事故下,水将射流进入铅铋环境,并与铅铋接触,引发剧烈传热、传质,同时铅铋与水蒸气密度相差甚大,难免会产生虚拟质量效应,因此本文考虑了虚拟质量力模型。本文所使用的虚拟质量力模拟参考了Fullmer等[19]开发的相关模型,并针对铅铋-水相互作用进行一定修正。

液态铅铋所受拽曳力模型为:

FVM,LBE=CVMαLρLαLBE,EαG·

(21)

式中:FVM,LBE为液态铅铋所受虚拟质量力矢量;CVM为相应系数;αL、αG分别为液体混合物、蒸汽混合物体积份额;ρL为液体混合物理论密度;αLBE,E、αH2O,E分别为液态铅铋、液态水在整个液体混合物中的相对体积份额;VG、VLBE、VH2O分别为蒸汽混合物、液态铅铋、液态水的速度矢量。

液态水所受拽曳力模型为:

FVM,H2O=CVMαLρLαH2O,EαG·

(22)

式中,FVM,H2O为液态水所受虚拟质量力矢量。

蒸汽混合物所受拽曳力模型为:

FVM,G=-CVMαLρLαG·

(23)

式中,FVM,G为蒸汽混合物所受虚拟质量力矢量。

上述虚拟质量力模型为矢量形式。另外由于Fullmer等提出的虚拟质量力模型针对的是气-水两相流,也是RELAP5/MOD3.3所使用的模型,而铅铋-水相互作用涉及多相流,针对该问题,本文对此模型进行了修正。首先上述模型中αL=αLBE+αH2O表示液相体积份额,考虑了液态铅铋和液态水两种液体,液相理论密度ρL按式(24)计算,且使用αLBE,E、αH2O,E进行了权重修正;其次对于液相加速度,由于铅铋-水相互作用条件下包含两种液体,在上述虚拟质量力模型中利用αLBE,E、αH2O,E进行了加权平均修正。

(24)

1.2.5交界面相变传热模型铅铋-水相互作用将引发剧烈的传热、传质,因此本文考虑了交界面相变模型,即质量守恒方程中的ΓC,H2O、ΓV,H2O等,由交界面传热确定,下面以液态水汽化为例进行详细阐述。

1) 确定交界面类型

2) 计算交界面净传热速率,判断相变类型

液态水-水蒸气交界面的净传热速率为:

(25)

3) 计算相变速率

如若交界面发生了汽化现象,则:

(26)

1.2.6状态方程模型状态方程(EOS)模型是封闭守恒方程的必要条件,同时对非平衡传热、传质迭代计算,压力迭代计算非常重要。铅铋-水相互作用涉及物质组分包括:液态铅铋、铅铋蒸气、液态水、水蒸气及不可凝结性气体,因此需要考虑以上5种物质组分的状态方程模型。

1) 蒸汽组分EOS模型

蒸汽组分包括铅铋蒸气、水蒸气,其状态方程参考Modified Redlich-Kwong(MRK)模型,如对于水蒸气而言,其基本形式如下:

(27)

(28)

式中:pG,H2O为水蒸气分压;RG,H2O为水蒸气气体常数;υG,H2O为水蒸气比体积;bG,h2o,1、bG,h2o,2、bG,h2o,3、bG,h2o,4分别为相应拟合系数;TH2O,CRT为水临界温度。

2) 不可凝结性气体

不可凝结性气体采用理想气体EOS方程模型。

3) 液体组分

液体组分包含液态铅铋、液态水两种,并考虑它们为可压缩流体,这里以水为例进行说明。液体组分温度按式(29)计算,式中TS,H2O为液体饱和温度,TS,H2O、∂TH2O/∂p均是液体比内能的单值函数,由拟合关系式计算。液体组分比体积υH2O按式(30)计算,式中υS,H2O为液体饱和比体积,υS,H2O、∂υH2O/∂p也是液体比内能的单值函数。

(29)

(30)

1.3 数值算法

描述铅铋堆SGTR事故下池内铅铋-水相互作用行为的模型复杂、行为演化急剧。如守恒方程中存在急剧的传热、传质行为,相界面演化行为,液态铅铋、液态水、蒸汽混合物3种速度场之间存在复杂的耦合关联性,速度场、压力场也存在紧密耦合关联性。因此,通过直接求解方案求解上述偏微分守恒方程十分困难,尤其是收敛性难以确保。考虑到上述问题,本文采用了时间步分割四步算法进行求解:1) 网格内源项计算;2) 网格间对流计算;3) 压力迭代;4) 时间步末更新。

1) 网格内源项计算

隐式联立求解质量守恒方程和能量守恒方程,仅考虑传热、传质源项。第1步的核心是基于EOS的非平衡态传热、传质迭代计算,迭代收敛后,更新全部密度和比内能。然后显式求解动量守恒方程,更新全部速度。

2) 网格间对流计算

先后分别显式求解质量、能量守恒方程,仅考虑对流项,第2次更新全部密度和比内能。然后,半隐式迭代求解动量守恒方程,考虑对流、压力、拽曳力及虚拟质量力项。其中,压力、对流项采用显式格式,拽曳力、虚拟质量力项采用隐式格式,迭代收敛后第2次更新全部速度,并获得速度对压力的偏导数。

3) 压力迭代

确定残差方程,共考虑5个:压力残差(EOS压力和网格压力之差)δp;液态铅铋、液态水、蒸汽混合物质量方程残差,即δm,LBE、δm,H2O、δm,G,仅考虑对流效应,采用半隐格式;蒸汽混合物能量方程残差δe,G,其中比内能对流采用显式格式,空泡份额对流项采用半隐格式。将上述残差方程对网格压力,液态铅铋、液态水、蒸汽混合物的宏观密度,蒸汽混合物温度,以及4个临近网格压力分别进行泰勒展开,即可获得压力迭代方程,并直接求解。最后结合EOS模型,开展循环迭代计算,直至上述5种残差收敛,即可更新压力、全部密度、蒸汽混合物温度及全部速度。

4) 时间步末更新

重复第2步,只是网格压力采用压力迭代最终值,对流项采用半隐格式,具体求解算法与第2步类似。第4步将更新全部变量,即获得时间步数值。

2 程序验证

2.1 日本JAEA水注铅铋实验

为验证所开发的程序模拟铅铋-水相互作用行为能力,首先选取了日本JAEA开展的水注铅铋实验进行验证,图2为实验装置示意图[5]。实验条件为:铅铋温度为778 K,压力为105Pa,水注射速度为4.7 m/s,注射水温度为298 K,注射孔口半径为3 mm,注射持续时间为500 ms。利用所开发程序对实验过程进行数值模拟,并与实验结果进行对比,如图3所示。由图3可看出数值模拟结果与实验结果符合较好。

图2 水注铅铋实验装置Fig.2 Experimental apparatus of water injection into LBE

2.2 意大利ENEA水注铅铋实验

日本JAEA水注铅铋实验主要聚焦于机理层面,而意大利ENEA为从铅铋堆SGTR事故工况角度研究铅铋-水相互作用行为,同时为验证SIMMER-Ⅲ程序,先后搭建了LIFUS5、LIFUS5/MOD2实验台架。本文选取其中两个实验用于程序验证:1) 基于LIFUS5实验台架的水注铅铋实验[11];2) 基于LIFUS5/MOD2实验台架的水注铅铋实验[12]。

1) 基于LIFUS5实验台架的水注铅铋实验

LIFUS5实验台架[11]含反应容器、膨胀容器、水池、排放池及液态金属储存池等。本文选取的实验条件为:液态铅铋温度为350 ℃,注射水压力为6×105Pa、温度为130 ℃,注射持续时间为10 s,反应容器内覆盖气体体积份额为5%。

本文建立了水注铅铋实验的二维轴对称数值仿真几何模型,如图4所示。通过模拟获得了反应容器内液态铅铋、覆盖气体的压力随时间的演化过程,并与实验结果进行对比,如图5所示。由图5可见,模拟结果与实验结果符合较好。另外还可发现在实验初期,反应容器内压力震荡比实验结果更剧烈,尤其是覆盖气体压力。

图4 LIFUS5实验台架水注铅铋实验数值仿真几何模型Fig.4 Geometrical model of simulating water injection into LBE experiment based on LIFUS5

图5 反应容器内压力的演化过程Fig.5 Variation process of pressure in reaction vessel

2) 基于LIFUS5/MOD2实验台架的水注铅铋实验

2015年,意大利ENEA对原LIFUS5实验台架进行了升级、改造,形成LIFUS5/MOD2[12]。该实验台架包含反应容器、加压水池、注射管线等。ENEA设置了4套实验方案,本文选取了第4套,实验条件为:铅铋温度为400 ℃,注射水压力为4 MPa、温度为240 ℃,反应容器内覆盖气体(氩气)体积份额约为23%,注射持续时间约为2.5 s,注射孔口直径为4 mm。

本文建立了该水注铅铋实验的二维轴对称数值模型,与图4类似,通过模拟获得了反应容器和加压水池内压力随时间的演化,并与实验结果[12]进行对比,如图6所示。由图6可看出,总体上无论是从定性角度,还是定量角度,模拟结果与实验结果均符合较好。就局部对比而言,在1.0～4.0 s期间,对于反应容器内压力,前后时间段内本文数值模拟压力均略高于实验值,中间时间段内本文数值模拟压力略低于实验值。