含空气蒸汽冷凝数值计算模型建立与模拟

宿吉强，孙中宁，张东洋

(1.哈尔滨工程大学 核安全与仿真技术国防重点学科实验室，黑龙江 哈尔滨 150001；

2.中国船舶重工集团公司 第703研究所，黑龙江 哈尔滨 150078)

新一代核反应堆系统广泛采用非能动安全壳冷却系统以保持安全壳的完整性。对于基于双层混凝土安全壳的非能动安全壳冷却系统，位于安全壳内的冷凝器是其最关键的设备之一，而安全壳内的不凝性气体(主要为空气)的运动与分布对冷凝器的影响至关重要。在蒸汽冷凝过程中，大量空气的存在会导致传热的严重恶化，所以必须对含空气的蒸汽冷凝开展研究[1]。

目前对于含不凝性气体蒸汽冷凝的研究主要有实验研究和理论分析两大类。实验研究得到的经验关联式具有较高的实用价值[2-5]，但经验关联式的适用条件苛刻，而具有广泛适用性的经典方法和理论尚有待形成。对于工程应用，通过数值模拟的方法获得评估数据是一种经济且快速的途径，同时数值模拟也可获得实验中较难得到的流场特性及动态特性的预测结果。

含不凝性气体蒸汽冷凝的数值模拟研究，主要利用CFD软件和计算机语言编程相结合的方式实现。目前主要的数值模拟方法有两类：一类是用基本的物理定律求解在冷凝表面的传热传质，在不凝性气体存在条件下建立一个两相边界层，从能量守恒和质量守恒原理出发提出一套数值计算模型，在数学方程的基础上解释不凝性气体影响蒸汽冷凝的机理[6]；一类是采用已发展的传热传质关系式的积分计算，将这些关系式转化为相应控制方程源项，应用于临近冷凝表面的单元网格内，从而实现冷凝过程的数值模拟[7]。第1类方法由于需要大量的编程及相当细致的表格，不适用于工程应用，而第2类方法可迅速直观地模拟冷凝过程中流场的变化情况，具有较高的工程应用价值，因而本文以此为基础对含空气的冷凝进行数值模拟。

为验证数值计算模型的准确性，并得到相应的流场分布结果，本文对Su等[8]和Cornet等[9]含空气蒸汽冷凝的实验进行数值仿真模拟，旨在为含空气蒸汽冷凝的传热传质研究奠定基础。

1 实验系统与实验方法

文献[8]的实验是在压力为0.2～0.6 MPa、空气质量分数为20%～80%的条件下进行的。实验模型示于图1。含空气蒸汽在长2 000 mm的竖直冷凝管的外表面进行冷却，混合气体存储在高3.55 m的碳钢圆柱罐体内，其直径为565 mm。蒸汽由罐体下侧入口射入，在入口与冷凝管之间存在两层均气孔板，冷凝管内通有冷却水。

2 计算模型

2.1 冷凝模型及理论假设

含空气蒸汽在竖直壁面冷凝的二维物理模型坐标系如图2所示。平行和垂直于冷凝壁面的坐标分别为y和x，对应的速度分量分别为v和u。主流气体速度为u，温度为T，空气质量分数为Wa，液膜层厚度为δf，空气积聚层厚度为δg，考虑浮生力，引入重力加速度g。

图1 实验模型简图

图2 坐标系及相关物理量

为了对冷凝过程进行有效的求解计算，本文对含空气蒸汽冷凝过程作了以下4点假设：

1) 冷凝过程中忽略液膜的影响，凝结液随着冷凝过程被移除，以研究不凝性气体在壁面附近的积聚及主流区的流场分布。

2) 混合气体中的蒸汽为其分压下饱和蒸汽，饱和蒸汽压力ps由饱和蒸汽温度Ts唯一确定，根据文献[10]，饱和蒸汽压力与饱和蒸汽温度的关系式为：

(1)

3) 冷凝仅发生在临近冷凝壁面的第1层网格内，气相边界层内饱和蒸汽受冷降温后，还有可能出现液滴。而参考文献[11]，与大量的壁面冷凝相比，气相边界层内液滴的存在对蒸汽的质量扩散仅产生极其微小的影响，因而对气相边界层中形成的液滴进行了忽略。

4) 将蒸汽空气混合流体处理成单相双组分流体，水蒸气和空气组成的双组分气体混合物假设为理想气体，推导其密度计算式为：

(2)

式中：ρmix为混合气体密度，kg/m3；R为摩尔气体常数，R=8.314 5 J/(mol·K)；Wv为水蒸气质量分数；Mv、Ma分别为水蒸气、空气摩尔质量，Mv=18 kg/mol，Ma=28.97 kg/mol。双组分混合物的质扩散系数由下式计算：

D=D0(T0/T)αp0/p

(3)

式中：T0、p0分别为标准状态下的温度、压力，T0=273.15 K，p0=101 325 Pa；D0为标准状态下质扩散系数，水蒸气和空气混合时其值为2.56×10-5m2/s；α取值范围为1.5～2.0，参考文献[12]取1.81。

2.2 含空气蒸汽冷凝的实现

图3 壁面单元示意图

基于前文的假设，可将含空气蒸汽冷凝过程简化为如图3所示的过程：在冷凝壁面第1层网格内，混合气体与壁面的对流换热导致流体温度降低，在冷凝壁面单元内比较当地蒸汽分压pv与壁面温度对应的饱和蒸汽压力ps的大小来判断是否发生冷凝，若pv>ps，则冷凝发生，产生凝液量m0，同时放出热量H0，凝液及潜热被从系统中移除；反之凝液量为零，仅存在对流换热。对流换热量由FLUENT计算得到，冷凝的相变过程会引起混合气相的质量、能量、动量及组分损失，这些损失可利用编译UDF(user defined function)在冷凝壁面的控制方程中添加源项来模拟，使用UDM(user defined memory)存储临近壁面各单元内的源项值，可用于计算某一时刻蒸汽释放的潜热。

模拟过程中，冷凝传热系数的获取采用TOSQON提出的理论[13]，将潜热与显热混合物通过气液界面向冷凝液传递的热量分为对流换热过程的显热和伴随凝结传质过程的潜热两部分。依据文献[14]定义m0的计算式为：

m0=CihoA(Tb-Tcw)/hlg

(4)

与之相对应的冷凝焓H0为：

H0=m0(cp,mixTb-cp,airTref)

(5)

式中：Ci为调整系数，W/(m2·K)；ho为实验得到的冷凝传热系数；A为近壁面网格单元面积，m2；hlg为饱和蒸汽的汽化潜热，J/kg；cp,mix、cp,air分别为混合气体和空气的比定压热容，J/(kg·K)；Tref为焓为0时的参考温度，Tref=273.15 K。

式(4)、(5)中除了Tcw，其余均为罐内主流的物理参数。选择Su等[12]提出的关联式作为ho的计算公式，调试获得的调整系数为0.98，得到m0的计算式为：

m0=0.98[10 189.3+90 416.4p-

(4 314.4+46 537p)lg(100Wa)]·

(Tb-Tcw)0.4/hlg

(6)

在模拟实验工况时，采用瞬态逼近稳态策略，选取基于压力的全隐式耦合算法。根据前文的假设，壁面冷凝只发生在近壁面单元。FLUENT并不能分辨出近壁面单元，因而需设定判断逻辑流程。为此，本文对求解近壁面单元冷凝的逻辑流程进行了设定：求解之前，通过壁面ID查找出近壁面单元的ID，并开辟存储空间储存这些ID，对流体域单元扫描计算控制方程时，判定该单元ID是否与存储空间ID相符合。如果符合，说明该单元是近壁面单元，进行冷凝计算；如果不符，说明该单元不是近壁面单元，不进行冷凝。需要注意，在每个时间步长的计算收敛后，需使用UDF保存壁面单元的压力和水蒸气质量分数及壁面温度，为下一时间步长是否发生冷凝提供依据。图4为UDF编译结构的流程图。

2.3 几何模型与处理方法

对二维几何模型进行网格划分，简化去除了冷却水的进、出口段，全部采用四边形结构化网格。均气孔板的小孔处利用尺寸函数加密，临近冷凝壁面的第1层网格高度为2 cm，通过网格无关性解的验证，确定二维模型的网格量为29 442。

图4 冷凝程序结构

在对模型进行计算时，对密度、动量、能量项和组分项均采用二阶迎风格式进行离散；选择二阶瞬态隐式关系式求解双组分气体混合过程的解；湍流模型选择RNGk-ε模型。

冷凝壁面采用无滑移壁面边界条件，给定壁面温度分布，冷凝壁温的差异会对传热产生较大影响[12]，因而模拟过程中采用实验的壁温参数。壁面温度通过拟合传热管轴向的外壁面温度的实验值得到。温度边界由壁温实验值的拟合关系式得到，并通过UDF的定义边界宏加载到冷凝壁面。除蒸汽入口外冷凝罐的其他壁面均按绝热处理。入口边界条件为质量流量入口，通过UDF定义入口质量流率等于UDM储存的每个单元凝液量的总和，用以模拟每一稳定工况的恒压过程。

3 数值模型验证及结果分析

3.1 模型验证

Su等[12]和TOSQON的参数范围[13]与安全壳事故工况下的参数相比，具有一定的代表性，因而本文在CFD中对该实验进行了模拟。

为验证本文方法的准确性，对TOSQON实验的含空气蒸汽冷凝的工况条件进行了数值模拟。TOSQON实验是在一高为4.8 m、直径为1.5 m的圆柱形压力罐内进行的，蒸汽及不凝性气体(空气、氦气)入射口位于距离罐底2.1 m处，入射口直径为41 mm，纯蒸汽入射时蒸汽流量为0.001 1 kg/s，入口蒸汽温度为126 ℃。冷凝壁面位于距离罐底2.391 m和4.391 m之间的区域，冷凝壁面温度为101.8 ℃，上、下两侧的非冷凝壁面温度分别为122.0 ℃和123.5 ℃，网格数为6 319。实验初始条件及模拟结果列于表1，初始条件误差较小。同时模拟结果显示，数值计算预测得到的温度、速度以及水蒸气含量同实验结果的误差分别为±1.89 K、±4.3%和±10%，体现出本工作计算方法在预测相似工况实验方面体具有较高的准确度。

表1 TOSQON初始实验条件及模拟结果

对Su等实验中3个压力的9组工况进行二维数值模拟，所选实验工况的主要物理参数列于表2。预测结果显示，冷凝传热系数最大相对误差为28.59%，主流温度及壁面过冷度的最大误差为3.75 ℃，空气质量分数最大相对误差为0.15%，进一步验证了本文所用模型的准确性。

表2 仿真计算结果与参考值对比

3.2 流场分析

图5为p=0.4 MPa、Wa=34.7%工况下的稳态模拟结果。在速度云图中，罐内的高速区形成于靠近冷凝罐和传热管的壁面附近，且冷凝罐下半部分空间的流速高于上半部分空间。在灌顶和罐底空间流场低速特点明显，说明在这两处流体较少参与罐内掺混，形成滞流区。流线图及局部矢量图显示冷凝罐内形成了多处涡流，从流线图可看出，传热管两侧的涡流形状修长，涡流中心位于传热管下端的左、右两侧，即在传热管两侧各自形成了自然对流。

a——速度云图；b——流线图；c——局部速度矢量

图6为p=0.4 MPa下不同空气质量分数时的温度分布云图。从图6可看出，随入射蒸汽量的减小，均气孔板以上空间的温度场不均性逐渐减弱，在传热管底端20 cm以内有低温流体的积聚。从图中的等值线可看出，在冷热流体交汇处形成了较大的温度梯度和明显的温度交界面，较高的温度使蒸汽和混合气体在混合过程中进行快速的热量交换，让蒸汽达到其当地分压下的饱和温度。高温气体沿冷凝罐内壁上升，同时与低温流体混合，在冷凝罐上部空间形成了均匀的温度场。在图6b中黑框所标出的位置，有两块温度高于周围流体的高温区，蒸汽沿远离入口侧的冷凝罐内壁面向上流动与罐内低温流体混合。结合图6a、c发现，受到下降冷流体在均汽孔板上部堆积的影响，蒸汽的上升混合过程是变化的，可能在靠近入口侧上升与罐内气体混合，也可能是在远离入口侧上升与罐内气体混合，或在两侧同时上升与罐内气体混合。

a——Wa=34.7%；b——Wa=57.2%；c——Wa=83.3%

4 结论

证明了利用实验结果实现含空气蒸汽冷凝数值模拟的可行性，该方法与其他模拟方法相比具有简单、可靠、快速等特点，具有一定的工程应用价值。利用本文模拟方法对TOSQON实验及Su等实验共10组工况进行了数值模拟，模拟结果同TOSQON实验值的相对误差在±10%以内，同Su等实验结果平均温度的最大误差为3.75 ℃，平均压力和空气质量分数的相对误差均在1%以内，计算结果与实际情况吻合较好，该数值计算方法对于相似工况以及相近物理模型内蒸汽在含空气的环境中冷凝过程的预测有指导意义。通过对数值模拟结果的分析，得到了速度场和温度场等参数细致的分布状态。

参考文献：

[1] ROSA J C. Review on condensation on the containment structures[J]. Progress in Nuclear Energy, 2009, 51(1): 32-66.

[2] UCHIDA H, OYAMA A, TOGO Y. Evaluation of post-incident cooling systems of light-water power reactors[C]∥Proceedings of International Conference on Peaceful Uses of Atomic Energy. [S.l.]: [s.n.], 1965: 93-102.

[3] TAGAMI T. Interim report on safety assessments and facilities establishment project[R]. Japan: Atomic Energy Research Agency, 1965.

[4] LIU H, TODREAS N E, DRISCOLL M J. An experimental investigation of a passive cooling unit for nuclear plant containment[J]. Nuclear Engineering and Design, 2000, 199(3): 243-255.

[5] DEHBI A A. Analytical and experimental investigation of the effects of non-condensable gases on steam condensation under turbulent natural convection conditions[D]. USA: Massachusetts Institute of Technology, 1990.

[6] ANDREANI M, PUTZ F, DURY T. Application of field codes to the analysis of gas mixing in large volumes[C]∥IAEA Technical Committee Meeting on the Implementation of Hydrogen Mitigation Techniques and Filtered Containment Venting. Koeln, Germany: [s.n.], 2001: 18-21.

[7] MARTIN J M, JIM M A, MARTIN F. Comparison of film condensation models in presence of non-condensable gases implemented in a CFD code[J]. Heat Mass Transfer, 2005, 41(11): 961-976.

[8] SU J Q, SUN Z N, FAN G M. Experimental study of the effect of non-condensable gases on steam condensation over a vertical tube external surface[J]. Nuclear Engineering and Design, 2013, 262(5): 201-208.

[9] CORNET P, MALET J, PORCHERON E. TOSQAN experimental results of the air-steam phase[C]∥OECD International Standard Problem on Containment Thermal-Hydraulics ISP-47. France: Institut de Radioprotection et de Sureté Nucléaire, Saclay, 2002.

[10] 刘志刚,刘咸定,赵冠春. 工质热物理性质计算程序编制及应用[M]. 北京：科学技术出版社,1992.

[11] 王朝阳,屠传经. 雾的形成对蒸汽和不凝性气体混合物的强迫对流膜状凝结换热的影响[J]. 高校化学工程学报,1989,3(1):75-84.

WANG Chaoyang, TU Chuanjing. The effect of fog formation on forced convective condensation of vapor noncondensable gas mixture[J]. Journal of Chemical Engineering of Chinese Universities, 1989, 3(1): 75-84(in Chinese).

[12] MASSMAN W J. A review of the molecular diffusivities of H2O, CO2, CH4, CO, O3, SO2, NH3, N2O, NO, and NO2in air and N2near STP[J]. Atmospheric Environment, 1998, 32(6): 1 111-1 127.

[13] COLBURN A P. Design of cooler condensers for mixtures of vapors with non-condensing uses[J]. Industry Engineering Chemistry, 1934, 26(11): 1 178-1 182.

[14] IVO K, MIROSLAV B. Modeling of containment atmosphere mixing and stratification experiment using a CFD approach[J]. Nuclear Engineering and Design, 2006, 236(14): 1 682-1 692.