GASFLOW程序液膜模型开发及初步验证

GASFLOW程序液膜模型开发及初步验证

王方年1，沈峰1，程旭1,2，黄兴冠2

(1.国核科学技术研究院，北京102209；2.上海交通大学，上海200240)

摘要：本文基于三维CFD安全壳程序GASFLOW开发了热构件壁面上的液膜覆盖与蒸发模型。通过选定AP1000大破口事故序列，采用耦合了液膜模型的GASFLOW程序分析了AP1000核电厂安全壳内温度压力响应及其非能动安全壳冷却系统(PCS)的性能，并与相同事故序列下WGOTHIC、MELCOR、CONTAIN等程序的计算结果进行比较。结果表明，耦合了液膜模型的GASFLOW程序可用于分析PCS的热工水力行为，其基本功能满足计算需要。

关键词：GASFLOW程序；液膜模型；AP1000；非能动安全壳冷却系统；安全壳响应

中图分类号：TL364.3 文献标志码：A

收稿日期：2014-07-16;修回日期：2015-01-22

基金项目：国家科技重大专项资助项目(2013ZX06004008)

作者简介：王方年(1985—)，男，江西九江人，工程师，硕士，核能科学与工程专业

doi：10.7538/yzk.2015.49.11.2044

Development and Preliminary Verification

of GASFLOW Code Coupled with Film Model

WANG Fang-nian1, SHEN Feng1, CHENG Xu1,2, HUANG Xing-guan2

(1.StateNuclearPowerResearchInstitute,Beijing102209,China;

2.ShanghaiJiaoTongUniversity,Shanghai200240,China)

Abstract:A model of heat structure wall surface film coverage and evaporation was developed based on 3 dimensional CFD containment code GASFLOW. The containment temperature and pressure response and the passive containment cooling system (PCS) performance of AP1000 during large break LOCA were analyzed by GASFLOW code coupled with film model. The calculation results were compared with the calculated results of other containment codes WGOTHIC, MELCOR, CONTAIN under the same accident scenario. The results show that the modified GASFLOW code coupled with film model is feasible to analyze the thermal-hydraulic behavior in PCS of PWR and the basic functions can meet the requirements for the calculation.

Key words：GASFLOW code; film model; AP1000; passive containment cooling system; containment response

GASFLOW程序是由美国LANL和德国FZK共同开发的三维计算流体力学程序[1-2]。它是一个最佳流场估计工具，可描述流场中以下三维流动的现象：气体的扩散、混合分布与分层，氢气燃烧和火焰扩散，不可凝气体的分布对本地凝结与蒸发的影响，气溶胶的夹带、传输与沉降等。

在AP1000核电厂安全壳设计中，由于采用了重力、压差和自然循环等非能动理念[3]，非能动安全壳冷却系统(PCS)与传统安全壳相比增添了许多新现象，尤其是钢制安全壳外壁面液膜的蒸发冷却现象。国内外研究人员曾使用多种程序(包括WGOTHIC、MELCOR、CONTAIN、COCOSYS等)对AP1000核电厂安全壳进行了分析[4-6]，但使用三维CFD程序分析的很少。基于AP1000核电厂安全壳的特殊性，GASFLOW能否适用于AP1000核电厂PCS分析需重新评估。GASFLOW3.2只能对所有构筑物表面输入初始液膜厚度，液膜不能持续不断的补充到钢制安全壳外壁面，因此液膜会随着时间推移而蒸干，无法实现安全壳的长期冷却。目前，最新版的GASFLOW3.5仍无法实现这一功能[1-2]。同时，GASFLOW用静态液膜厚度(液膜厚度为常数)来实现蒸发，无法表示液膜厚度随流动方向慢慢变薄的过程，也不能区分干湿区域，无法模拟液膜在钢制安全壳外壁面覆盖率的变化。

本文开发热构件壁面上的液膜覆盖与蒸发模型，并耦合到GASFLOW3.2中，对大破口事故下AP1000核电厂安全壳响应和PCS排热性能进行分析。

1液膜模型开发

液膜模型开发主要分两部分。首先是对GASFOW定义的钢制安全壳壁面网格进行重构，改变钢制安全壳表面的液膜计算顺序，使其能模拟液膜流动的走向；然后通过耦合液膜覆盖与蒸发模型对具有一定计算顺序的每一个网格上的液膜进行计算处理，求解液膜的质量守恒方程、能量守恒方程，计算出液膜覆盖率、排热能力等。

1.1GASFLOW网格重构

GASFLOW中热构件壁面网格按照式(1)从小至大的顺序进行编号，记为master[1]：

i+mboff(iblk)

(1)

其中：i、j、k为定义网格时的空间坐标编号；mboff(iblk)为描述建模对象的编号。故在GASFLOW中，各墙体壁面计算顺序并不是模拟AP1000液膜流动所需的先后计算顺序，因此，本文对GASFLOW安全壳壁面网格进行重构，以一定的方式存储，然后以一定的顺序提供程序计算，简单地实现液膜从里至外、自上而下流动的模拟。

将GASFLOW中每个壁面的网格进行了重新定义，通过一4维数组master=(i，j，k，windex)将所有壁面网格存储起来，然后按从里至外、自上而下进行循环计算，该数组包含每个墙体具体位置及其网格面朝向等信息。其中：i、j、k定义了空间坐标x、y、z3个方向的网格编号，即定义了具体位置；windex记录着每个网格面的朝向，本文用1、2、3来分别代表x、y、z方向，用-1、-2、-3来分别代表-x、-y、-z方向。

1.2液膜覆盖与蒸发模型

在完成GASFLOW壁面网格重组基础上，通过开发液膜覆盖与蒸发模型，计算不同位置、不同时刻的水膜厚度、水膜覆盖率；同时也能从液膜厚度的角度模拟溪流变窄过程及蒸干现象；最终较为准确地实现在安全壳壁面上由内到外的热量传递。

本文是基于Nusselt液膜理论分析液膜流动情况的，其基本假设是：液膜沿竖直方向匀速流动，根据重力与摩擦力等力的平衡求解速度。计算时认为1个网格上的液膜速度是均匀的，且是液膜厚度的函数；沿流动方向上的下一个网格，因为蒸发作用厚度产生变化，其速度也跟随变化，但网格上的液膜速度是同一值。

本文的液膜覆盖与蒸发模型只是针对钢制安全壳外壁面液膜开发的，实际情况中液膜也同时存在于钢制安全壳内壁面。钢制安全壳内壁面的液膜计算采用GASFLOW中原有模型，即静态液膜模型，认为液膜厚度是受限的，液膜通过冷凝达到一定厚度后将保持不变，继续凝结的部分将被移除出内壁面。

1) 平均液膜速度与液膜厚度

利用Nusselt理论[7]求解竖直壁面层流降膜的液膜厚度及液膜速度。动量方程为：

(2)

边界条件为：

(3)

定义竖直壁面单位宽度上的质量流量为：

竖直壁面上的液膜厚度为：

(4)

管外壁降膜流动的液膜厚度为：

其中：δpl为竖直平板的液膜厚度；εR为管外壁面液膜厚度与管道外径的比值。安全壳属于大直径的管外降膜流动，其外壁面液膜厚度(μm级)远小于安全壳直径(40 m或以上)，即εR≈0，可认为δ=δpl，因此式(4)仍适用于安全壳外壁面液膜流动。式(4)在推导时液膜流态属于层流，但实验证明[8]也符合湍流计算，甚至在相当一部分区域比层流计算更准确。Nusslet理论关系式在层流及波动层流的范围内计算出的液膜厚度均大于其他经验公式，液膜越厚可导致覆盖率的减小，因此本文采用Nusslet理论关系式计算是保守的。

2) 临界液膜厚度

西屋公司的LST、SST、Flat Plate实验[4]表明，当液膜流速大于某一临界值时液膜宽度保持不变，当液膜流速降至该值后流速保持不变，液膜宽度减小直至完全蒸发。由于润湿面积越大，蒸发流量越多，因此较大的流速限值使得液膜宽度更早减小，是保守的。本文基于文献[4]的液膜蒸发实验结果，取具有包络性的保守值作为临界液膜线质量流量参考值，再代入式(4)可得到临界液膜厚度。

3) 液膜覆盖率及蒸发流量的计算

定义m为液膜质量流量(即未蒸发的液膜流量)，W为湿周(液膜湿润宽度)，Z为壁面网格的高度(或长度)，则有[4]：

m=WΓ

(5)

液膜流动主要分为两个过程，即液膜宽度不变的过程(W不变，液膜厚度随蒸发变薄)和单位长度上的质量流量不变的过程(Γ不变，液膜已经到达临界厚度，并保持不变)。

总的液膜蒸发流量为液膜初始化总流量与未蒸发流量之差，即：

(6)

其中：mevp为液膜蒸发流量；minitial为液膜的初始注入流量；mn为第n个网格的液膜流量。

每一网格上液膜蒸发流量是采用传热传质类比的形式计算得到的[2]。当液膜向下流动并蒸发，液膜厚度减薄至低于临界液膜厚度后，钢制安全壳外壁面网格上的液膜将出现干、湿两部分区域，在程序模型内认为该网格的液膜宽度变窄，液膜高度不变，原理如图1所示。

图1　液膜覆盖率变化原理图 Fig.1　Schematic of variation for film coverage rate

由于液膜模型能计算出每一个网格的液膜宽度，因此液膜覆盖率Rcoverage可表示为：

(7)

其中：Wn、Zn分别为第n个网格上液膜的宽度与高度(或长度)；W0、Z0分别为网格的宽度与高度初始化值。每个网格上的液膜蒸发流量为：

(8)

(9)

因液膜蒸发带走的能量可表示为：

(10)

其中，hfg为液膜汽化潜热。

气体与热构件表面对流换热量为[10]：

(11)

其中：A为网格面积；Tbulk,gas为大空间混合气体的温度；Tsurf,structure为热构件表面的温度。

辐射换热量qradiation计算采用的是壁面-气流的简单辐射换热模型(壁面假设是黑体)：

(12)

其中：σ为Stefan-Boltzmann常数；εbulk,gas为气体发射率；αbulk,gas为气体吸收率。

2程序初步验证

本文以分析AP1000大破口事故下安全壳响应为例，初步验证耦合了液膜模型的GASFLOW程序的适用性与可信度。

2.1AP1000建模

1) AP1000 GASFLOW网格划分及几何建模

根据AP1000核电厂安全壳的实际布置情况进行了GASFLOW建模[3-4]。AP1000核电厂安全壳体积庞大且壳内结构复杂，若要建立一个与实际结构完全一致的模型则需采用很精细的网格，这将使计算代价变得无法接受。为寻求计算代价与计算精确度的平衡，在保持自由容积、传热面积的参数相对准确的情况下对安全壳构筑物进行了简化。径向网格由墙体和设备位置进行定位，轴向网格由各层地板和设备位置进行定位，尽可能保证网格界面与墙体、地板、设备表面重合以节省网格数目；圆周向采用均匀网格划分。本文安全壳模型采用的是圆柱坐标系下结构化正交网格划分方法。表1列出AP1000的网格划分方式及其数目，平均网格体积为1.18 m3。图2示出AP1000核电厂安全壳的GASFLOW模型。AP1000 GASLOW模型中安全壳内自由容积(不包括外部自然循环空间)为59 200 m3，钢制安全壳的高度为65.6 m，直径为39.6 m。

表1　AP1000 GASFLOW网格划分

2) AP1000大破口事故质能释放源项

AP1000大破口事故下，质能释放位置是在操作平台之下的蒸汽发生器隔间，质能释放源项参见文献[5]，包含：1) 从压力容器侧喷放出来的两相流体；2) 从蒸汽发生器侧流出的高温蒸汽；3) 自动卸压系统ADS-4释放在蒸汽发生器隔间的蒸汽。

3) 其他初始边界条件

安全壳内空气及热构件的初始温度为48.89 ℃，空气初始压力为0.108 MPa，除内置换料水箱隔间外，相对湿度为0。外部环境大气初始温度为46.11 ℃，压力为0.101 MPa，相对湿度为22%。钢制安全壳外壁面液膜蒸发的驱动力为界面与大气内水蒸气的浓度差，湿度用于计算大气内水蒸气的浓度。钢制安全壳与混凝土屏蔽构筑物之间的环形通道自然循环计算，本文采用定压边界条件。

图2　AP1000安全壳的模型 Fig.2　Geometric model of AP1000 containment

2.2大破口事故分析

为验证修改后GASFLOW的计算结果是否正确，与文献[5，11-12]的计算结果进行对比分析。表2列出AP1000大破口事故下安全壳内压力、温度的对比。本文计算结果与WGOTHIC、MELCOR、CONTAIN等集总参数程序的计算结果相近，但略有差别。AP1000核电厂安全壳的设计压力为0.508 MPa，设计温度为148.89 ℃[3]，本文的计算结果再次证明AP1000当前PCS的设计满足大破口事故下的安全需求。

表2　大破口事故下安全壳压力和温度的结果对比

图3示出相同边界条件下AP1000大破口事故的安全壳内压力的对比。压力呈现双峰值形式，前一个峰值是由于压力容器侧短时间内(1～20 s)喷放的大量高温、高压两相流体造成的；在20～100 s间，安全壳内压力变得较为平坦，其原因是这一阶段压力容器侧喷放量减少至较低的水平；在100～1 700 s期间，壳内压力持续上升，其原因一方面是压力容器侧释放的流体焓值上升，另一方面是SG侧开始释放蒸汽。337 s后PCS投入，压力曲线此时出现拐点。在1 700 s之后，压力因SG侧蒸汽释放结束而缓慢下降。可见PCS的作用一方面是为安全壳的长期排热提供最终热阱，另一方面降低安全壳的峰值压力。

图3　不同程序计算的安全壳内压力对比 Fig.3　Pressure comparison inside containment by different codes

从图3各程序计算结果对比可见，压力变化趋势基本一致，集总参数程序间的结果较为接近， GASFLOW计算结果较集总参数程序计算结果差别较大。这是由于各程序本身及建模过程存在差异，如GASFLOW采用传热传质类比的形式，集总参数程序一般采用冷凝经验关系式。

图4示出不同程序计算的安全壳内温度的对比。集总参数程序一般情况下将安全壳内穹顶当作一个结块或几个结块处理，取结块的平均值；GASFLOW数据取自如图2所示的A、B、C3点温度的平均值。安全壳穹顶的平均温度均在事故前期质能高速喷放阶段(20 s内)达到峰值，高温流体会因喷放速度的影响迅速冲到安全壳穹顶，并随着热构件对高温蒸汽凝结的作用，使得穹顶气流温度达到峰值后又下降，而后因PCS的投入带走了热量使得温度缓慢下降。从长期结果看，4种程序的计算结果处于同一水平且具有相同的变小趋势。

图4　安全壳内穹顶处温度对比 Fig.4　Dome temperature comparison inside containment

图5　液膜覆盖率的变化 Fig.5　Variation of coverage rate for film

2.3PCS排热性能

图5、6分别示出钢制安全壳外壁面液膜覆盖率的变化及不同高度处液膜厚度的变化。图6中，a～g的位置如图2所示。由图5、6可见：因PCS喷淋尚未启动，液膜覆盖率在0～337 s时为零；PCS启动后，液膜覆盖率的变化与液膜厚度的变化相对应，且大部分监测点(a、b、c、d、e)的变化趋势保持一致；有些区域喷放前期没有液膜覆盖而后出现了液膜覆盖，并处于临界液膜厚度或蒸干状态(f、g处)。可见，本文程序关于液膜覆盖与蒸发模型的耦合具有一定合理性。

图6　钢制安全壳外壁面液膜厚度的变化 Fig.6　Variation of film thickness on outside wall surface of steel containment

图7示出钢制安全壳外壁面热流密度的变化。由图7可见，液膜蒸发平均热流密度约在11 kW/m2，PCS环形通道内对流换热、辐射换热的平均热流密度分别在0.5 kW/m2、0.2 kW/m2左右。因此，液膜蒸发是安全壳内热量导出的主要方式。

图7　钢制安全壳外壁面平均热流密度的变化 Fig.7　Variation of average heat flux on outside wall surface of steel containment

图8示出钢制安全壳外壁面不同监测点处温度的变化。钢制安全壳外壁面液膜温度超过100 ℃的时间不算长，超过100 ℃的区域集中在穹顶部分，竖直壁面未超过限值。本文程序暂时不能处理液膜温度超过100 ℃的情况，这是未来需改进之处。

图8　钢制安全壳外壁面温度的变化 Fig.8　Variation of outside wall surface temperature for steel containment

3结论

1) GASFLOW程序的二次开发满足预设的计算功能，能够计算液膜覆盖率、钢制安全壳外壁面液膜厚度的变化，能够模拟分析AP1000核电厂PCS的排热性能。

2) 对比分析相同事故序列下不同程序的计算结果，初步验证程序修改后的计算结果是可信的。

3) 通过与WGOTHIC等程序计算结果的对比，再次表明目前AP1000核电厂PCS的设计在大破口事故下能够保持安全壳的完整性，满足设计要求。

参考文献：

[1]SPORE J W, TRAVIS J R, ROYL P, et al. GASFLOW: A computational fluid dynamics code for gases aerosols, and combustion, Volume 2: User’s manual, LA-13357-MS, FZKA-5994[R]. USA: LANL; Germany: FZK, 1998.

[2]TRAVIS J R, SPORE J W, ROYL P, et al. GASFLOW: A computational fluid dynamics code for gases aerosols, and combustion, Volume 1: Theory and computational model, LA-13357-

MS, FZKA-5994[R]. USA: LANL; Germany: FZK, 1998.

[3]林诚格. 非能动安全先进核电厂AP1000[M]. 北京：原子能出版社，2008.

[4]WOODCOCK J, ANDREYCHEK T, CONWAY L, et al. WGOTHIC application to AP600 and AP1000, Rev. 1, WCAP-15862[R]. America: Westinghouse, 2004.

[5]TILLS J, NOTAFRANCESCO A. Application of the MELCOR code to design basis PWR large dry containment analysis, SAND2009-2858[R]. America: SANDIA, 2009

[6]PHILIPP B, HANS-JOSEF A. Simulation of AP1000’s passive containment cooling with the German containment code system COCOSYS[C]∥Nuclear Energy for New Europe. Bovec, Slovenia: [s. n.], 2011.

[7]SINKUNAS S, GYLYS J. Analysis of laminar liquid film flowing down a vertical surface[C]∥Fourth International Conference on CFD in the Oil and Gas. Norway: [s. n.], 2005.

[8]YU Y Q, WEI S J, YANG Y H, et al. Experimental study of water film falling and spreading on a large vertical plate[J]. Progress in Nuclear Energy, 2012, 54: 22-28.

[9]BIRD R B, STEWART W E, LIGHTFOOT E N. Transport Phenomena[M]. New York: John Wiley and Sons, 1960.

[10]ROHSENOW W M, CHOI H. Heat, mass, and momentum transfer[M]. New Jersey: Prentice-Hall, 1961.

[11]OFSTUN R. AP1000 LOCA containment pressure analysis, Rev.2, APP-SSAR-GSC-523[R]. America: Westinghouse, 2003.

[12]Westinghouse. Shield building air inlet flow area design change, Rev. B, APP-GW-GEE-1307[R]. America: Westinghouse, 2010.