PISAA软件中含不凝性气体的蒸汽冷凝模型研究与验证

李贺, 杨小明, 张宁娜, 马如冰, 元一单

(中国核电工程有限公司 中核核电安全严重事故研究重点实验室, 北京 100840)

核电厂事故期间一回路会向安全壳内释放出高温高压的气体,使壳内温度压力骤然升高,为保证安全壳在极限事故工况下的完整性,必须对其内部进行降压冷却。安全壳壁面上的冷凝传热是导出内部热量的重要手段之一[1-2],然而由于混合气体中常含有不可凝气体(如空气等),一定程度上阻碍了蒸汽在壁面的冷凝,极大增加了蒸汽在壁面冷凝的传热热阻,从而降低了蒸汽冷凝传热系数[3-4],可能导致热量不能及时带出安全壳,进一步威胁到反应堆的安全。

因此,含不可凝气体的蒸汽冷凝过程是一个非常复杂的过程,涉及到蒸汽与冷凝液体之间的质量、能量和动量的交换传递问题,蒸汽与不可凝气体之间的动量、能量交换和扩散问题。国内外许多学者对此现象进行了大量理论和实验研究[5-11],得到许多具有较高实用价值的理论分析结果和实验经验关系式,但由于相关实验关系式适用范围有限,仍没有具有广泛适用性的经典方法和理论。工程上通过数值模拟的方法既能有效获得实验中的流场和动态特性,又是一种经济高效的获得评估数据的方法[12]。因而合理恰当地使用计算模型、如实地模拟相关现象、提高冷凝传热计算的准确度对于核电厂事故条件下安全壳热工水力分析程序来说至关重要。

当前核电厂严重事故分析中主要采用美国开发的MAAP和MELCOR,以及欧盟开发的ASTEC等一体化软件,这些软件的核心内容完成于20世纪欧美核电主导期,有些软件架构和数值方法已经过时,有些物理模型(如冷凝)使用的经验关系式也存在对一些特殊工况的过渡外延,计算的科学准确性有待进一步验证,亟待进行模型更新。中核集团中国核电工程有限公司在综合考虑国际上主流的严重事故一体化软件的程序架构、物理模型以及国际上最新的严重事故研究成果的基础上,开发了适用于“华龙一号”堆型的严重事故一体化分析软件(program integrated for severe accident analysis,PISAA),该软件采用集总参数算法,将连续的流体空间分割成一块一块的控制体,并通过流道相互连接起来,将壁面模化成热构件,同时综合考虑了现有堆型的创新设计,使用更加广泛、更加通用的计算模型,更加现代化的高效计算方法,能够快速对压水堆核电机组的严重事故主要现象做出准确判断和推演。

本文针对严重事故一体化分析软件PISAA中的冷凝计算模型,通过构建具体计算工况,对PISAA软件中含不凝性气体的蒸汽冷凝模型进行了对比验证,结合公开的Wisconsin冷凝实验数据[13]和ISP-47中的TOSQAN实验数据[14]对模型计算结果做了进一步确认,并对计算结果进行了分析。

1 冷凝传热模型

1.1 传热机理

当固体壁面温度低于与其相接触的含有水蒸气的气相的露点温度时,冷凝现象将会发生。图1为含不可凝气体蒸汽在竖直壁面的冷凝模型,由于蒸汽冷凝在壁面上形成的液膜沿重力方向流动,不可凝气体含量增加,且随着液膜的增厚,其流动状态由层流转变成湍流,当液膜厚度较薄时,对流效果不明显,液膜热阻也可以忽略处理。同时沿y轴负方向,蒸汽浓度减小,不可凝气体浓度增大,其分压大于混合主流气体中的不可凝气体分压,从而形成反向扩散驱动力使不可凝气体向主流扩散;而主流混合气体中,蒸汽分压高于气液界面处的分压,为蒸汽穿过不可凝气体层到达壁面发生冷凝提供了驱动力,蒸汽向冷凝表面的扩散与不可凝气体向主流气体方向的扩散维持着总压不变的动态平衡过程。因此对于含有不可凝气体的凝结过程,不可凝气体在气液交界面处的逐步汇聚,增加了传热传质的阻力,使得冷凝传热系数相比纯蒸汽冷凝时大幅度降低。

图1 含不可凝气体蒸汽在竖直壁面的冷凝模型[3]Fig.1 Condensation model of vapor containing non-condensable gas on vertical wall[3]

理论研究主要有2类研究方法：第1类是求解边界内的连续方程、动量方程、能量方程及气体组分方程得到相应的解析解,考虑了两相流、层流和湍流结构及液膜波动的影响,但是由于高度复杂和高额花费,使得这些模型在核安全分析中并不适用;第2类是基于扩散理论求解扩散边界层内的传热传质方程。此类方法可针对不同工况灵活采用不同条件下的计算模型,通过对冷凝传热各过程的独立计算,最终获得相关的热工水力参数。此类方法可以针对不同工况灵活采用不同条件下的计算模型,因而目前大多分析程序多数采用此种方法计算相变。本文介绍的PISAA软件中对冷凝计算的数值模拟计算也采用了第2类方法。

1.2 传热传质比拟模型方法

此方法主要是基于扩散理论,通过使用相似原理来求解扩散主导的传质速率。对于含有不可凝气体的相变计算,目前使用较为广泛的是 Kreith[15]模型和 Collier[16]模型。这2种模型理论基础基本一致,都是基于菲克定律,考虑了含有不可凝气体作用下发生的蒸汽的质扩散和法向主流蒸汽分量2个部分构成的冷凝传质过程,并沿边界层积分,从而得到总物质的量通量密度。

根据Colburn等[17]提出的经典传热传质比拟方法,对于含不可凝气体的蒸汽传热分析,蒸汽/非凝性气体边界层传递的热量主要包括蒸汽和非凝性气体的显热以及蒸汽冷凝所释放的潜热2个部分,其等效热阻可以表示成图2所示。

图2 等效热阻示意Fig.2 The equivalent thermal resistance of the wall condensation process

根据能量守恒定律：

q/htot=q/hfilm+q/(hconv+hcond)

(1)

式中：hfilm为液膜传热系数;hconv为对流传热系数;hcond为冷凝换热系数。将式(1)展开,可得：

hcond(Tb-Ti)+hconv(Tb-Ti)=hfilm(Ti-Tw)

(2)

式中：Tb、Ti、Tw分别为主流气体温度、相界面处气体温度以及壁面温度。

由于hfilm、hconv、hcond的计算均需要液膜表面温度Ti,因而为了求得总换热系数htot,需要通过假设液膜表面温度初值,通过关系式(2)得到新的Ti,通过不断迭代,直到新旧相界面温度Ti之差达到可以接受的偏差,从而得到最终液膜表面温度Ti。得到Ti后即可求得总换热系数htot。

1.3 液膜传热系数hfilm计算

通常,Nusselt理论解以及基于Nusselt理论的修正方程广泛应用于计算液膜传热系数,然而对于Nusselt理论,其在假设液膜为层流流动,当竖直表面过长时,随着高度的降低,液膜充分发展,液膜内流动由层流逐渐转变为波状层流,进而变为湍流,单纯采用根据Nusselt理论解会产生一定偏差,在Reδ≤30时采用文献[18]：

(3)

30

(4)

在Reδ>1 800时,采用Labuntsov建议的关系式[18]：

(5)

(6)

计算出Reδ后,最终计算出液膜传热系数hfilm[18]：

(7)

1.4 对流换热系数hconv计算

按对流换热相似准则数关联式,计算不考虑法向传质影响下的对流换热系数hconv：

hconv=Nu·k/l

(8)

式中Nu根据对流形式的不同参考外掠平板的湍流自然对流和强制对流关系式[14]：

NuFC=0.037Re0.8Pr1/3

(9)

NuNC=0.13(Gr·Pr)1/3

(10)

对于混合对流,采用Churchill关系式[19-21]：

(11)

2种流动方向相同取正号,相反取负号,对于不确定流向时,采用保守估计取最小值。

1.5 冷凝传热系数hcond计算

对于凝结过程,有：

qcond=hcond(Tg-Ti)

(12)

(13)

因而：

(14)

(15)

式中：hconv为对流换热系数;λ和T为扩散边界层气体热导率和温度;Dsteam为蒸汽的扩散系数;Msteam为蒸汽摩尔质量;P为气体总压;Psteam,i、Pstm,b分别为蒸汽在相界面和主流处的分压;Pnon,avg为气相不可凝气体平均分压;Pr、Sc为普朗特数及施密特数。计算中假定相界面处水蒸气为饱和蒸汽,水蒸气在壁面发生膜状凝结,冷凝后液膜均匀的附着在壁面上,同时根据冷凝质量及冷凝水的密度,结合壁面传热面积,可表示出冷凝液膜的厚度[22-23]。

2 Wisconsin实验验证

针对热构件表面冷凝传热模型的验证,本文选取了与Wisconsin常压部分冷凝实验结果进行对比[5]。Wisconsin实验主要研究含不可凝气体的大空间平板壁面冷凝传热现象,实验考虑大空间内竖直和水平壁面冷凝2个部分,容器右上方顶部和侧壁分别排列6块总长为0.914 4 m、厚度0.304 8 m的铝制冷凝板,不可凝气体分数为0.31～0.8,压力为105Pa,容器内部为定温定压的湿空气,湿空气在铝制冷凝板壁面发生冷凝传热[13]。实验时,冷凝板外侧被持续地过冷水冷却,从而使冷凝板壁面温度保持恒定,而大空间容器内通过从容器底部喷入高温蒸汽实现并保持初始不同工况下的定温定压状态。通过设定不同的初始热工参数,用HFM(heat flux meters)和CEB(coolant energy balance)2种方式测量各实验工况下的冷凝传热系数。表1给出了各实验工况下必要的初始参数。

表1 冷凝传热初始工况Table 1 Initial condition of condensing heat transfer

Wisconsin实验中在测量对流传热系数时使用了2种不同的方式,但是并未给出2种测量结果的不确定度,因而为提升实验对比结果的可靠性,对比时采用2种测量方式获得结果的平均值。

建模时,冷凝板壁面一侧对流传热,另一侧保持定温,大空间内的湿空气在壁面处发生凝结传热。针对不同初始工况分别按照PISAA软件输入卡要求制卡,图3展示了竖直和水平冷凝板冷凝传热系数的实验测量值与程序计算值对比情况。

图3 实验测量与PISAA计算冷凝传热系数对比Fig.3 Comparison of condensation heat transfer coefficient by experimental measurement and PISAA code calculation

从图3中可以看出,初始各工况下的PISAA计算冷凝传热系数与实际实验测量值基本一致,其中竖直冷凝板算例中,工况9～11的计算值与实验测量值偏差略大,最大在29.147 7%,但仍在工程可接受的偏差范围内。而水平冷凝传热系数计算值与测量值偏差小得多,平均偏差在2.548 1%左右,整体与实验测量值符合较好。对比结果表明,PISAA程序中使用的冷凝模型能够覆盖不同摩尔分数、不可凝气体情况下多元化壁面温度和主流温度的冷凝传热工况,对冷凝板方向不同而引起的结果差异也能够准确模拟。因此,PISAA软件中使用的冷凝传热计算模型科学有效且合理正确。

3 ISP-47 TOSQAN实验验证

TOSQAN容器由不锈钢制圆柱形腔室和一个地坑构成。内部总容积(含地坑)为7.0 m3,高4.0 m,圆柱最大直径1.5 m,壁面温度通过2个导热油回路被控制为2个恒定的温度,圆柱形腔室分为顶部、中间区域和底部3个部分,其中顶部、底部和地坑称为“热壁面”或“不冷凝壁面”,温度相同;中间区域称为“冷壁面”,为冷凝区,高度为2 m,传热面积9.4 m2,可控制在不同的温度值,地坑的高度为0.87 m,最大直径0.68 m,蒸汽、空气和氦气沿管线注入腔室内部,喷嘴直径0.041 m,高度为2.1 m,居中置于内部,几何形状信息见图4所示[14]。

图4 TOSQAN台架几何参数[14]Fig.4 TOSQAN geometry parameter[14]

TOSQAN实验的实验工况如图5所示,整个过程分为阶段A和阶段B 2个部分,共分为4个稳态,其中稳态1、稳态2和稳态3无不可凝气体的注入,属于阶段A;稳态4是在注入不可凝气体氦气后达到,属于阶段B。

图5 ISP-47 TOSQAN实验工况Fig.5 ISP-47 TOSQAN overview of test sequence

首先,在实验开始前1 d保持壁面恒温,实验开始后,首先进行第1次加压,将蒸汽注入干燥的热空气中,注入流量为1.4～1.14 g/s,随时间渐变,容器内部压力逐渐升高且达到恒定;然后向蒸汽注入管线中以恒定流量添加示踪粒子(用于激光测速的固体气溶胶颗粒),持续10 min,系统达到第1个稳态;第1个稳态持续一段时间后,系统进行第2次升压,蒸汽注入流量增加至12.27 g/s,系统经过短暂过渡后达到并保持在第2个稳态;然后通过降低蒸汽注入流量来降低系统的压力到达第3个稳态,与第1个稳态压力相同,仅初始条件不同;第3个稳态持续一段时间后,再次向系统中添加示踪粒子,以补偿在壁面发生沉降的气溶胶,随后沿蒸汽注入管线向系统中注入氦气,注入时间持续600 s,注入流量为1.03 g/s,此过程蒸汽流量从1.11 g/s随时间线性变化降低至0.89 g/s,而后系统逐步达到第4个稳态。

实验初始条件中系统绝对总压为105Pa,上部非冷凝壁面平均温度为122.0 ℃,下部非冷凝壁面平均温度为123.5 ℃,冷凝壁面平均温度为101.3 ℃,平均气体温度为115.4 ℃,空气体积分数为100%,无氦气及蒸汽。

图6为ISP-47 TOSQAN实验壳内压力测量值。从图6的时间坐标轴可以看出,TOSQAN ISP-47实验历时约11 h。压力的升高伴随着冷凝传质流量的提高,直到后者升高至与入射蒸汽流量相等,系统达到平衡。稳态1和稳态2为了建立不同入射蒸汽质量流量下的平衡态,从而验证不同入射流量下的程序计算能力。稳态3是为了验证程序计算中不同初始条件的影响,与稳态1唯一的不同是：稳态1初始TOSQAN台架内并没有蒸汽分布,而对于稳态3,TOSQAN台架内已经存在了蒸汽。稳态4是为了验证对于含轻质不可凝气体对流动、质量分布以及膜态冷凝现象的数值模拟。

图6 ISP-47 TOSQAN实验壳内压力测量值Fig.6 ISP-47 TOSQAN volume pressure by measured

本文根据TOSQAN实验的几何参数及实验工况,结合PISAA软件计算输入卡编制要求,制作了此实验的计算输入卡。计算采用集总参数模型,将整个系统划定成一个控制体,将容器内的上、中、下冷凝板及地坑分别设定成热构件,边界条件根据各时间段气体流量及温度的不同分别设置,并假定前一时间步长计算过程中产生的冷凝水对后续阶段的冷凝无影响。

图7展示了PISAA软件计算的ISP-47 TOSQAN实验过程中系统压力的变化情况,从图中可以看出,系统压力随着蒸汽的注入逐渐升高,之后过热蒸汽在热构件壁面出现冷凝,冷凝水与注入的热蒸汽质量上达到动态平衡,系统达到第1个稳态,后随着蒸汽注入量升高,系统压力在较大幅度增加后再次达到动态平衡状态即稳态2,蒸汽注入量再次降低时系统压力也在逐步降低,第3个稳态后同步注入蒸汽与不可凝气体,可以看到,不可凝气体的注入使得系统再次达到稳态的耗时要稍微长一些,也说明了不可凝气体的存在影响了蒸汽的正常冷凝。

图7 PISAA计算值与实验测量值对比Fig.7 Comparison between PISAA calculation and experimental data

整体来看计算结果与实验测量值吻合较好,能够准确模拟实验各阶段的压力变化情况,对于系统压力变化的瞬态区间亦能够如实展现,整个过程中,稳态2计算值与实验测量值偏差最大,在5%左右,其余阶段基本无偏差。

4 结论

1)针对事故工况下核电厂安全壳内热构件冷凝传热问题,PISAA软件采用的基于扩散边界层内传热传质相似理论(HMTA)的Kreith模型,模型科学合理,计算方法得当准确,能够准确模拟含有不可凝气体的蒸汽冷凝现象。

2)PISAA软件对Wisconsin冷凝传热实验的计算验证中,水平冷凝传热的传热系数结果平均偏差在2.548 1%,竖直冷凝板的计算偏差较大,但也仍在工程可接受的范围内,造成偏差较大的原因可能是竖直壁面中冷凝液膜厚度的简化方式与实际情况不太一致,后续可对此做进一步的优化和改进。

3)根据ISP47-TOSQAN的初始条件设定,PISAA软件能够基本还原实验过程的压力变化,且压力对比值的最大偏差出现在稳态2阶段,约5%左右,其余阶段基本无偏差,后续将进一步细化控制体节点,对冷凝水的质量及传热量等参数做详细分析。相关计算结果说明软件在处理含不可凝气体的冷凝传热过程中方法正确合理,能够满足安全壳内的热工水力计算需求。