基于Realizable k-ε 湍流模型的氨气泄漏数值模拟研究*

孙恩吉，李红果，王　敏

(1. 中国安全生产科学研究院，北京 100012；2. 国家安全生产监督管理总局信息研究院，北京 100012)

0　引言

随着我国工业化生产和存储的快速发展，现代化工艺的液氨制冷冷库不断增多，冷库无论从规模上、设备上和自动化程度上都有了显著提高。然而，我国仍存在很多老旧冷库，这些冷库由于设备陈旧、安全管理落后导致生产安全事故时有发生，安全生产形势依旧严峻。2012年10月22日，湖北省洪湖市德炎水产公司发生氨气泄漏事故[1]，造成400多人氨气中毒抢救。2009年8月5日，内蒙古赤峰制药集团氨水配制车间液氨卸车过程中，发生氨气泄漏[2]，致使202人接受治疗。液氨制冷站制冷系统所需液氨储存在储罐及输液管道内，液氨输运管道、阀门或储罐容器壁等处意外破损、爆裂都会导致液氨大量泄漏，液氨泄漏后会迅速蒸发并在站内迅速扩散[3]，对站内及附近人员的人身安全将造成致命的危害。

采用计算流体动力学软件来数值模拟泄漏扩散的过程，研究气象条件、地表情况、泄漏源等对气体扩散的影响是气体泄漏扩散研究的重要方法之一。熊立春等[4]引入时间叠加的高斯液氨泄漏扩散模拟及人员疏散；J. Labovsky等[5]对Fladis试验过程进行数值模拟，使得CFD动态模拟得到成功应用；Galeev等[6-7]模拟了氨气扩散的过程，研究了连续释放过程中风速和障碍物对扩散分布范围的影响；S.M. Tauseef等[8]计算了在障碍物情况下重气的扩散方式；郑茂辉等[9-10]着重分析了建筑物的扰动和不同风速对氨气浓度的影响；王志鹏[11]采用了高斯模型对液氨储罐泄漏事故进行了研究，提出了泄漏后的快速处置方法；于加收[12]对不同风向条件下氨的泄漏扩散规律进行了研究；吴其荣等[13]模拟了火电厂区的氨泄漏及其扩散状况。以往数值模拟更多的关注液氨储罐事故与自然风影响下的区域扩散规律，针对室内不同的泄漏点、不同的排风条件下氨气的运移规律的研究相对较少。

因此，提出了基于Realizablek-ε的氨气泄漏限元数值模拟分析方法，计算了液氨泄漏质量，分析了增加排风口、不同液氨泄漏口及不同液氨泄漏量的氨气扩散规律及浓度变化情况，对于设置应急处置装备、采取合理措施具有指导意义。

1　模型的建立

1.1　氨气泄漏物理模型

氨气泄漏物理模型包括液氨储罐，排风口A，B，进风口，空气压缩机3台。液氨制冷车间尺寸为20 m×9 m。液氨储罐为高5 m，直径2.5 m的圆形罐体。压缩机尺寸为1.5 m×1 m。进、排风口直径尺寸均为1 m，分别位于左侧的上下方及右侧上方；通风口的尺寸为1 m，模型具体布置如图1所示。

网格划分是流体数值模拟的重要基础，在模型离散化的过程中，模型被划分成若干离散量元素。运用六边形和非结构化的网格划分，将模拟区域划分成3部分不同网格密度的计算区域，其中，泄漏口上方和中间泄漏区域采用最密集的四边形网格，两侧区域网格为较疏的四边形网格，网格包括大约500 000六边形网格元素和400 000四边形网格元素。以下关于液氨泄漏蒸发后的氨气扩散过程模拟，选取其上部排风口为压力出口,下方通风进口为速度入口边界条件。

1.2　氨气泄漏量的计算

在常压下，液氨在空气中会极速膨胀，大量气化，并扩散到较大的空间范围。液氨制冷车间正常运转时，液氨储罐泄漏发生在液相区域，在常温、常压下将以气态形式存在，假设泄漏路径较短，不会形成汽化核心而使部分液氨在泄漏路径上汽化形成闪蒸两相流，同时，假定液氨泄漏时全部被蒸发，而不考虑形成液池，其泄漏速率用伯努利方程计算：

(1)

式中:Q0为液氨泄漏速率，kg/s；Cd为泄漏系数；A为泄漏口横截面积，m2；Pt为贮氨器内压，Pa；P0为大气压，Pa；h为泄漏点距液面高度，m；圆形泄漏口，其Cd取值为0.65，泄漏口直径为0.01 m；240 K时，液氨密度ρt为681 kg/m3；储罐内压Pt为950 kPa；大气压力P0为101 kPa；h取值为0.4 m；g取值为9.8 m/s2；液氨泄漏质量初始流量为1.645 kg/s。

泄漏的液氨会发生闪蒸，泄漏时直接蒸发的液体所占百分比F可按下式计算：

(2)

式中：CP为液体的定压比热，J/kg·K；T为泄漏前液体的温度，K；T0为液体在常压下的沸点，K；H为液体的汽化热，J/kg。根据经验，当F>0.2时，一般不会形成液池；当F<0.2时，F与蒸发液体之比，有线性关系，即当F=0时，没有液氨被蒸发；当F=0.1时，有50%的液体被蒸发。

以下数值模拟主要针对制冷车间内液氨储罐短时间(10 s，30 s和60 s)泄漏的氨气运移规律，根据制冷车间内部的实际布置和工况，可作如下基本假设：

1)假设液氨储罐环境温度与储罐初始温度相同。

2)液氨储罐泄漏口处的氨气泄漏速度恒定。

3)液氨储罐泄漏过程中氨气作为理想气体考虑。

4)泄漏过程不发生任何相变反应和化学反应。

5)进风口处的风速恒定。

6)泄漏过程中，温度恒定且与外界无热量交换。

1.3　Realizable k-ε模型

Realizablek-ε模型由Launder and Spalding提出后，出现了RNGk-ε模型及Realizablek-ε模型，被广泛应用于工业流场及热场的交换模拟中。Realizablek-ε模型为湍流粘性增加了1个限制条件公式，并为耗散率增加了传输方程，对旋转流动、流动分离及复杂二次流都有更好的实现，对平板和圆柱射流发散比率有更精准的预测。利用Realizablek-ε模型在对氨气泄漏过程中的射流现象的数值分析更能体现出优越性，其表达式如下：

(3)

(4)

2　模拟结果分析

2.1　增加排风口的浓度分布

根据《呼吸防护用品的选择、使用与维护》(GB/T18664-2002)中的规定，空气中氨气立即威胁生命或健康的浓度(IDLH)为360 mg/m3。选择此浓度为临界点，换算成体积摩尔浓度为2.12×10-5kmol/m3。模拟将排风口B设置为压力出口，其他条件不变，增加排风口A时，计算泄漏时间为10 s，30 s，60 s时的氨气浓度变化，结果如图2所示。

图2　氨泄漏10 s, 30 s, 60 s时氨气IDLH区域(增加排风口A)Fig.2　Ammonia leak 10 s, 30 s, and 60 s when ammonia IDLH area (increased exhaust port A)

由图2分析可以得出，增设排风口A，氨气泄漏10 s后，氨气由泄漏口垂直排入泄漏区，因其比空气轻，上升至一定高度并开始在顶层飘逸。由于排风口A侧下方为进风口，氨气仍偏向排风口B侧进行扩散，并扩散至泄漏区的顶部，在其到达顶部后，受到顶部阻力的影响，开始向下并沿四周空间扩散。在氨气扩散和上升过程中，其高浓度区域不连续，并出现“扰流”现象。氨气泄漏30 s后，泄漏氨气向排风口B侧扩散的更快速，B处氨气浓度也相对较高，泄漏至顶部的氨气沿顶部扩散的范围和区域也有所增加，并且在顶部方向向排风口A侧扩展。氨气泄漏60 s后，泄漏氨气气流开始向排风口A侧扩散，进风口对扩散的影响也逐渐加大。氨气在到达B处顶部后，开始沿顶部及右侧墙壁向下进行扩散。

2.2　不同泄漏口的浓度分布

将泄漏口位置设置在储罐右上部，边界条件设置为速度入口，泄漏口泄漏量为1.645 kg/s，将原泄漏口设置为固壁边界，其他条件不变，模拟泄漏口位置不同、泄漏时间为10 s，30 s，60 s时的氨气浓度变化情况，结果如图3所示。

图3　氨泄漏10 s, 30 s, 60 s时氨气IDLH区域(改变泄漏口位置)Fig.3　Ammonia leak 10 s, 30 s, and 60 s ammonia IDLH area (changing leak port position)

改变泄漏位置，氨气泄漏10 s后，泄漏氨气气流向排风口B侧方向扩散，由于泄漏时间较短、泄漏速度大，氨气并未向上扩散。氨气泄漏30 s后，泄漏气流开始向排风口B侧扩散，受到第2台压缩机的阻挡，氨气向厂房顶部扩散，储罐与压缩机间的氨气浓度开始增大，氨气在到达排风口B处浓度也开始逐步增加。氨气泄漏60 s后，进风口的风速对泄漏扩散的影响增大。由于储罐与压缩机间无法进行气体扩散，其间的氨气浓度泄漏接近泄漏口处的氨气浓度，氨气气流向排风口B侧扩散趋势明显。

2.3　不同泄漏量的浓度分布

将泄漏量由1.645 kg/s分别增加至1.81 kg/s，1.99 kg/s，其他条件不变时，模拟泄漏区泄漏量变化，以及泄漏时间为10 s，30 s，60 s时的氨气浓度变化情况。泄漏时间为10 s，泄漏量为1.81 kg/s，1.99 kg/s时氨气浓度变化对比如图4所示。

图4　氨泄漏10 s时氨气IDLH区域Fig.4　Ammonia leak 10 s ammonia IDLH area in different Leakage(1.81 kg/s upper panel, 1.99 kg/s lower panel)

泄漏量增加、泄漏时间10 s时，泄漏区氨气最高摩尔浓度由0.397 kmol/m3增加至0.399 kmol/m3；泄漏量为1.99 kg/s时，氨气浓度更容易上升至顶部，同时也主要集中在排风口B侧。泄漏时间为30 s，泄漏量为1.81 kg/s，1.99 kg/s时氨气浓度变化对比如图5所示。

图5　氨泄漏30 s时氨气IDLH区域Fig.5　Ammonia leak 30 s ammonia IDLH area in different Leakage(1.81 kg/s upper panel, 1.99 kg/s lower panel)

泄漏区氨气最高摩尔浓度由0.420 kmol/m3增加至0.426 kmol/m3，泄漏口竖直方向氨气高浓度区域较多，排风口B处的氨气浓度升高。泄漏时间为60 s，泄漏量为1.81 kg/s，1.99 kg/s时氨气浓度变化对比如图6所示。泄漏氨气扩散到达上部顶板后，沿顶板向四周运动。

图6　氨泄漏60 s时氨气IDLH区域Fig.6　Ammonia leak 60 s ammonia IDLH area in different Leakage(1.81 kg/s upper panel, 1.99 kg/s lower panel)

3　结论

1)通过数值模拟分析发现，在增加排风口的情况下，氨气泄漏60 s后，与进风口同侧的排风口氨气浓度增加幅度小于1%，氨气主要朝对流方向的排风口飘逸。在不同泄漏口模拟条件下，由于氨气气流会受到厂房内不同设备设施的阻挡，60 s后进风口的风速对氨气泄漏扩散影响增大。不同氨气泄漏量的情况下，通过1或2个湍流模型可产生更为合理的数值模拟结果。

2)数值模拟结果的重要性是所选择的模型需要考虑在多相局部稳定性的影响和介质以及涡流扩散的非各向同性效应。氨气的分子量小于空气的平均分子量，对流排风口位置对于泄漏氨气浓度影响较大，增加对流端的排风口有利于氨气气云的快速扩散和稀释。氨气泄漏量的增加使得泄漏口垂直方向氨气浓度增加显著，同时，空气流动速度增加导致氨气扩散趋于稳定，能有效减小泄漏源造成的灾害危险范围。

3)基于Realizablek-ε湍流模型的氨气泄漏有限元数值模拟可以有效的分析和预测氨气泄漏后运移规律，该方法可为液氨制冷企业应急救援预案的制定、氨泄漏事故预防及应急处置提供参考。

[1] 张兵.冷库建筑防火设计分析——以宝源丰禽业有限公司特别重大火灾事故为例[J].中国安全生产科学技术, 2013, 9(10):157-167.

ZHANG Bing. Analysis on fire safety design of cold storage building based on the lesson of fire disaster in BAO YUAN FENG Poultry Co. Ltd [J]. Journal of Safety Science and Technology, 2013, 9(10):157-167.

[2] 庄学强.大型液化天然气储罐泄漏扩散数值模拟[D] .武汉：武汉理工大学，2012.

[3] 赵明,陈求稳.城市重大危险源区域风险评价方法研究[J].中国安全生产科学技术, 2014, 10(9):158-164.

ZHAO Ming, CHEN Qiuwen. Study on method of regional risk assessment for urban major hazard [J]. Journal of Safety Science and Technology, 2014, 10(9):158-164.

[4] 熊立春,陈建宏,石东平.引入时间叠加的高斯液氨泄漏扩散模拟及人员疏散研究[J].中国安全生产科学技术, 2015, 11(11):76-82.

XIONG Lichun, CHEN Jianhong, SHI Dongping. Simulation on liquid ammonia leakage by Gaussian model introducing temporal superposition and personnel evacuation [J]. Journal of Safety Science and Technology, 2015, 11(11):76-82.

[5] J. Labovsky, L. Jelemensky. Verification of CFD pollution dispersion modeling based on experimental data [J]. Journal of Loss Prevention in the Process Industries, 2011, 24(2): 166-177.

[6] A. D. Galeev, E. V. Starovoytova, S.I. Ponikarov. Numerical simulation of the consequences of liquefied ammonia instantaneous release using FLUENT software [J]. Process Safety and Environmental Protection, 2013, 91:191-201.

[7] A. D. Galeev, A. A. Salin, S.I. Ponikarov. Consequence analysis of aqueous ammonia spill using computational fluid dynamics [J]. Journal of Loss Prevention in the Process Industries, 2013, 26: 628-638.

[8] S.M. Tauseef, D. Rashtchian. S.A. Abbasi. CFD-based simulation of dense gas dispersion in presence of obstacles [J]. Journal of Loss Prevention in the Process Industries, 2011, 24(4): 371-376.

[9] 郑茂辉，金敏，许建明.城市建筑群环境有毒有害气体扩散数值模拟[J].同济大学学报自然科学版，2013, 41(1):48-52.

ZHENG Maohui, JIN Min, XU Jianming. Numerical simulation of hazardous gas dispersion around buildings in urban environment [J].Journal of Tongji University (Natural Science) 2013, 41(1):48-52.

[10] 郑茂辉，郭月容.城市居住区毒气扩散数值模拟[J].环境保护科学，2012, 38(1):1-4.

ZHENG Maohui, GUO Yuerong. Numerical simulation of toxic gas dispersion in urban residential areas [J].Environmental Protection Science, 2012, 38(1):1-4.

[11] 王志鹏.氨气泄漏数值模拟及应急响应措施研究[D]. 北京；中国地质大学，2013.

[12] 于加收.开放空间液氨泄漏扩散规律及人员疏散的研究[D].北京：中国地质大学(北京)，2014.

[13] 吴其荣，杜云贵，喻江涛，等. 火电厂氨泄漏扩散特性模拟研究[J].环境科学与工程.2014, 37(2): 155-159.

WU Qirong，DU Yungui，YU Jiangtao, et al. Simulation of characteristics of ammonia leak diffusion in thermal power plant [J].Environmental Science & Technology, 2014, 37(2):155-159.