周书航, 边浩志, 吴桐雨, 李文涛, 丁铭
(1.哈尔滨工程大学 核动力装置性能与设备黑龙江省重点实验室,黑龙江 哈尔滨 150001; 2.哈尔滨工程大学 核科学与技术学院,黑龙江 哈尔滨 150001)
当反应堆发生破口事故时,会向安全壳内喷放大量的高温高压蒸汽,致使安全壳内压力和温度急剧上升乃至威胁其完整性[1-2]。为了确保安全壳边界完整性,非能动安全壳热量导出系统普遍应用在第三代核电厂中,其主要利用流体密度差产生驱动力,通过蒸汽冷凝的方式持续带走安全壳内热量[3-5]。然而在蒸汽冷凝过程中,安全壳内不可避免地会存在不凝性气体[6-9]。
文献[10]研究表明,不凝性气体的存在会极大削弱蒸汽冷凝能力,即使只存在1%的空气,冷凝传热系数也会下降60%。随后国内外学者对含不凝性气体蒸汽冷凝现象开展了大量的实验研究[11-13]以及数值分析[14-17],其中,数值分析中普遍采用的模型是扩散边界层模型。该模型可以在较广的参数范围预测蒸汽冷凝局部流动和传热现象,但在高压、高蒸汽浓度下的适用性仍然存在问题。
针对这一问题,Bird[18]、Dehbi[19]、Bian[20]等分别提出了修正因子进一步提高了扩散边界层模型的准确性和适用性。但目前国内外在对含不凝性气体蒸汽冷凝过程进行数值分析时,由于不凝性气体的存在,几乎忽略了液膜热阻。然而,在高压、高蒸汽浓度下液膜热阻是否可忽略,以及液膜将如何影响蒸汽冷凝特性尚未揭示。因此,本文采用CFD 计算软件STAR-CCM+(17.02)开展含空气蒸汽冷凝和层流液膜耦合模型研究,评估液膜效应如何影响蒸汽冷凝特性以及在不同参数范围内的变化机制。
含空气蒸汽冷凝-层流液膜耦合模型的建立需要针对多组分气体侧、液膜侧以及气-液交界面分别建立控制方程描述其冷凝的完整过程。
气体侧控制方程主要描述多组分气体的流动与传热过程,其中控制方程主要包括质量守恒方程、动量守恒方程、能量守恒方程和组分输运方程。
质量守恒方程:
(1)
动量守恒方程:
ρgfg+Sg,pw
(2)
能量守恒方程:
(3)
组分输运方程:
(4)
式中:ρ为流体密度,kg/m3;t为时间,s;w为速度矢量,m/s;Sm为质量源项,kg/(m3·s);P表示表面力,N/m2;τ为剪切力,m2/s;f为体积力,N/m3;Spw为动量源项,N/m3;E为流动流体所具有的能量,J;ω为质量分数;keff为流体的等效传热系数,W/(m·K);Sh为能量源项,J/(m3·s);D为扩散系数,m2/s;下标i表示组分,下标g表示气体。
此外,气体侧流速较高时,其流态会由层流过渡到湍流态,需在此基础上添加湍流模型。本文采用可实现的k-ε两层湍流模型。该模型适用于各种流动的过程,相比于标准k-ε模型有更好的预测精度,具有适用性广、精度高的优点。
层流液膜控制方程主要描述液膜的层流流动与传热过程,其中控制方程主要包括质量守恒方程、动量守恒方程、能量守恒方程和组分输运方程。
质量守恒方程:
(5)
动量守恒方程:
(6)
能量守恒方程:
(7)
组分输运方程:
(8)
式中:δ是液膜厚度,m;下标l表示液体。
气-液交界面控制方程主要描述蒸汽冷凝传热过程,其冷凝传热过程采用基于扩散边界层理论的冷凝模型,该模型可以很好地预测局部流动和传热现象。当蒸汽在气-液交界面冷凝时,气相和液相之间将发生质量、动量和能量交换。
质量源项:
(9)
式中Δ为近壁面第1层网格厚度,m。其中,混合气体的扩散系数D计算公式为[21]:
(10)
式中:下标0表示标准状态。
动量源项:
Spw=Sm·w
(11)
能量源项:
Sh=Smhv
(12)
为评估基于扩散边界层建立的含空气蒸汽冷凝-层流液膜耦合模型,需开展模型验证工作。COAST实验[6]的参数范围(压力:0.2~1.6 MPa,空气质量分数:0.16~0.71,过冷度:45~117 ℃)较广,相比于安全壳事故工况,其具有一定的代表性。为在较广的参数范围内评估模型,本文针对COAST实验工况进行了选取并建立了数值模型。该模型为直径1.5 m、高度2.5 m的竖直壳体,壳体中心设置有长度1 m,直径19 mm的单管,如图1所示。冷凝段设置为恒壁温壁面,入口采用质量流量入口。
图1 几何与网格模型示意Fig.1 Geometry and mesh model
为确保该模型网格条件的计算精度,需开展网格无关性验证工作,由于本文所采用几何模型与作者前序研究所采用的模型相同,因此网格条件的选取(主流网格尺寸m= 0.04 m 和Y+= 1,网格数量为30万)与前序研究保持一致[20]。
通过对COAST实验工况数值模拟发现,如图2所示,随着压力、蒸汽质量分数的提高,扩散边界层模型预测的蒸汽冷凝传热系数与实验值偏差逐渐变大,尤其在高压、高蒸汽质量分数条件下,偏差可达到50%以上。主要原因考虑是扩散边界层模型无法预测冷凝率增强引发的抽吸效应以及液膜流动对气侧影响的液膜效应,因此需对模型进行修正。
图2 扩散边界层模型验证Fig.2 Validation of diffusion boundary layer model
由于扩散边界层模型本身的缺陷,Bird认为当蒸汽在换热面发生冷凝后,换热面附近和主流区域之间会产生较大的气体浓度梯度,从而增强气体之间的对流扩散以强化换热,这种现象称为抽吸效应[19-20]。
为了在模型中考虑这一额外对流/扩散传质项,本文通过自定义场函数的形式在冷凝模型中添加了Bird修正因子B和θB,并将扩散系数DBird添加到模型中:
(13)
(14)
(15)
DBird=D0·θB
(16)
此外,Dehbi认为Bird修正因子会高估其冷凝率,因此基于Bird修正因子提出了新的修正因子θD,并与Dehbi实验进行了验证,关系式为:
θD=0.346 64+0.495 24θB+0.165 02θB2
(17)
本文在模型中分别添加Bird和Dehbi修正因子与实验值进行对比,如图3所示。在模型中添加Bird和Dehbi修正因子后,其预测的传热系数偏差虽相较于扩散边界层模型来说已大幅减小,但整体偏差大于20%。
图3 各种修正因子比较Fig.3 Comparisons among various correction factors
其主要原因考虑是Bird和Dehbi修正因子均不是在自然对流条件下提出,且都是基于忽略液膜的假设,而为使模型适用于自然对流条件下分析且考虑液膜效应,需对其进一步修正。
由于该模型考虑了液膜的影响,当处在高压、高蒸汽浓度工况下时,蒸汽冷凝率会显著提升,液膜厚度也随之增加,导致更大的液膜热阻以及液膜轴向流速,如图4所示。
图4 液膜厚度及流速变化规律Fig.4 Variation of liquid film thickness and velocity
然而,液膜热阻增大会导致换热系数降低,液膜轴向流速增大会提升冷凝换热系数,因此液膜效应对冷凝换热既有强化作用,也有抑制作用。为更好地预测自然对流条件下的液膜效应,本文基于COAST实验提出新的修正因子θZ:
(18)
为评估新修正因子θZ,本文基于COAST实验、Dehbi实验工况[22]进行了计算对比验证,如图5、6所示。结果表明:与Bird和Dehbi修正因子相比,新修正因子θZ与实验结果吻合较好,偏差普遍在10%以内。
图5 修正后模型验证Fig.5 Validation of modified model
图6 耦合模型预测值与实验值对比Fig.6 Comparison between predicted results of the coupled model and experimental results
由于目前CFD研究中通常在含不凝性气体蒸汽冷凝条件下忽略液膜热阻,本文为评估忽略液膜产生的偏差以及液膜效应的作用机制,在P=1.6 MPa,ΔT=20 ℃工况,对不同空气质量分数下忽略液膜的单管外含空气蒸汽冷凝进行了数值计算,如图7所示。结果表明:忽略液膜后,当空气质量分数高于0.5时,偏差小于10%,而当空气质量分数小于0.5,偏差将迅速上升,在空气质量分数ωa=0.1时偏差达到20%以上。由图8也可以看出,随着空气质量分数的减小,总热阻和气侧热阻呈现下降的趋势,而液膜热阻呈现上升的趋势,并且在空气质量分数小于0.5时,液膜的热阻占比会快速增加,最高可达60%以上。基于上述结果得出,当处在高压、高蒸汽浓度工况时,液膜热阻的影响是不可忽略的,有必要在模型中考虑液膜效应的影响。
图7 忽略液膜偏差分布Fig.7 Deviation distribution under ignoring liquid film
图8 不同空气质量分数下热阻分布Fig.8 Thermal resistance distribution under different air concentration
此外,本文评估了液膜效应在不同工况参数下的作用机制,如图9所示。结果表明:忽略液膜模型将明显低估实际传热过程中的换热面气膜轴向流动速度,并且随着压力、蒸汽质量分数的增大,其模型低估导致的偏差将越来越大。基于以上分析,液膜对换热既有强化作用也有抑制作用,强化作用主要是液膜会带动近壁面气膜起到轴向加速的作用,从而增强冷凝换热。抑制作用主要是液膜本身存在热阻导致。这一结果进一步证实了液膜的影响是不可忽略的,尤其是在高压、高蒸汽浓度条件下。
图9 有液膜与无液膜下气体流速对比Fig.9 Comparison of gas velocity with and without liquid film
在高压高蒸汽浓度条件下,冷凝量会大幅提升,从而使得冷凝液膜急剧变厚,除了液膜热阻变得不可忽略外,其液膜的流态可能不限于层流状态。为获取实验中液膜流态区域(层流区域、湍流区域),本文针对实验参数范围进行了数值计算,提取对应工况下的液膜最大雷诺数,如图10所示。结果表明:COAST实验中压力P≤ 1.3 MPa的实验工况,其冷凝液膜流态均处在层流区域,模型可较好地预测其局部流动与换热特性。
图10 不同工况下最大雷诺数Fig.10 Maximum Reynolds number under various conditions
当处于压力P=1.6 MPa、空气质量分数ωa=0.16工况,其冷凝液膜最大雷诺数达到3 000,表征其冷凝液膜在部分高度已处于湍流区域。为了获取液膜湍流区域长度,图11示出了压力P=1.6 MPa、空气质量分数ωa=0.16时不同过冷度下的局部雷诺数分布。结果表明:在P=1.6 MPa、ωa=0.16工况下,其液膜层流区域长度约为0.6~0.7 m,处在换热管头部;湍流区域长度约为0.3~0.4 m,处在换热管尾部。
图11 P = 1.6 MPa、ωa=0.16下局部雷诺数分布Fig.11 Local Reynolds number distribution at P=1.6 MPa, ωa=0.16
1)提出的含空气蒸汽冷凝-层流液膜耦合模型,可在广的参数范围内较好地预测自然对流条件下的含空气蒸汽冷凝过程,其偏差普遍在10%以内。
2)冷凝液膜对换热既有强化作用,也有抑制作用。抑制作用是由于其本身存在热阻,强化作用主要是液膜对近壁面的气膜具有轴向加速作用。当空气质量分数低于0.5时,忽略液膜会产生较大偏差,最大可达到20%以上。
3)总结得到了COAST实验的液膜流态区域。在高压、高蒸汽浓度条件下,管外冷凝液膜会达到湍流态。当处在压力P≤1.3 MPa的实验工况,其冷凝液膜流态均处在层流区域,当处于压力P=1.6 MPa工况,其尾部液膜已处于湍流区域。