时婷婷，杨福明，常雪伦，孙 晨，于文涛，鲁仰辉

(国家电投集团科学技术研究院有限公司，北京 102209)

0 引言

氢能作为21世纪的绿色能源，得到世界各国的普遍关注。随着氢能的快速发展，氢安全问题正受到越来越多的关注。相比于其他能源形式，氢具有一系列独特属性。氢分子是自然界中最小的分子，在高压系统中很容易通过管道接头、阀门等位置发生泄漏；氢对很多材料具有氢脆效应，容易导致临氢环境下的承压设备失效；氢气密度低、扩散快，在露天场合不易形成可燃混合物，但是在受限空间内则很容易在局部发生聚集；同时，氢气具有可燃范围宽、燃烧热值高、燃烧速度快、爆炸能量大的特点，聚集的可燃气云一旦被点燃可能会发生燃爆事故，其后果往往非常严重。综上，在开展氢能及其应用技术的基础研究过程中，需明确氢安全事故演化行为及后果，探索针对氢安全事故所展开的应急处理措施，为相关法律法规及标准提供可靠的数据支撑，扎实推动氢能技术发展和规模化应用。

泄漏往往是氢安全事故的始发事件，也是后续燃爆行为分析的输入条件，因此氢气泄漏扩散行为及演化规律是氢安全研究的基础性内容。针对受限空间内氢气泄漏扩散的规律，国内外研究学者利用实验和数值模拟进行了相应的研究。Denisenko 等[1]和Manins[2]通过实验研究总结了氢气在不同形状、大小和边界条件的受限空间内的扩散演变机制。实验结果显示：在相同泄漏流速的情况下，氢气在受限空间内的扩散行为可分为2 种，当氢气的初始泄漏速度小于150 m/s 时，氢气在受限空间内的扩散方式为填充模式；而当氢气的初始泄漏速度大于等于150 m/s 时，氢气在受限空间内的扩散方式变为渐淡模式，从安全角度来讲，渐淡模式比填充模式更危险。

de Stefano 等[3]通过实验分析了受限空间内泄漏位置及泄漏速度对氢气扩散的影响。实验结果显示：氢气在受限空间内的扩散行为主要取决于泄漏速度而不是泄漏的初始位置，泄漏的初始位置对氢气扩散的影响随着泄漏速度的增加而减小。当泄漏速度较小时，氢气首先聚集在受限空间的底部，然后再上升至天花板，最后下降到其他位置，直至均质化阶段。

氦气的安全性高，可操作性强，且氦气与氢气同为低密度气体，氦气扩散实验可作为与在受限空间内使用氢气相关的安全研究的一部分，在泄漏扩散这一物理过程中，充分了解氦气这种低密度气体在空气中的分散机制至关重要；并且由于氦气的安全性高，用来模拟氢气在受限空间内的扩散行为的可操作性相对较强。Cariteau 等[4]通过实验分析了在1 m×1 m×1 m 的受限空间内3 种不同通风口形状下氦气的扩散行为。实验结果显示：氦气的扩散行为受通风口形状的影响较大；自然通风会显著改变低密度气体的扩散行为。

根据上述研究结论，可在氢能基础设施中设置便于通风的开口。Pitts 等[5-6]、Blais 等[7]通过13 组1∶1 全尺寸实验，研究了车库内氢气泄漏扩散及燃烧的发展过程。实验结果为氢燃料电池车辆相关标准和规范的开发提供了氢相关安全数据，以减轻氢燃料电池的火灾隐患。随后，Hajji 等[8]通过Fluent 软件对文献[5-7]研究的车库内氢气泄漏实验进行了模拟，对比验证了泄漏位置、释放速率、泄漏持续时间和氢气浓度之间的关系，模拟结果与实测结果吻合良好。

准确的数值结果可为氢能设施的传感器布设、结构优化和通风装置设计等提供参考，为此，基于OpenFOAM[9]、EPPL 算法[10]和多组分库Mutation++[11]，本项目研发团队研发了鲁棒、高效、精确的全速域多组分浮力流求解器——MPPBuoyantEpplFoam 求解器。本文采用MPPBuoyantEpplFoam 求解器对氢气在受限空间内(即室内环境)的泄漏扩散行为进行了数值模拟，对氢气泄漏后的流动扩散规律开展了研究，并对其多组分扩散功能的准确性进行了验证。

1 模型描述

1.1 实验装置

本文以de Stefano 等[3]的实验装置作为模拟对象。实验场景为某核电厂内易发生氢气泄漏的隔间，为了便于调试氢气流速，实验装置为该隔间1/15 缩比模型；实验装置的形状为立方体，长×宽×高的尺寸为0.47 m×0.33 m×0.20 m；泄漏口位于装置顶部。为了探究不同位置氢气浓度随时间的变化关系，实验装置内设有16个传感器，分别固定在A、B、C、D 4 根支架上，每根支架上有4个传感器。实验装置的尺寸及布置的传感器位置示意图如图1所示，传感器的具体位置坐标如表1所示。

表1 传感器的具体位置坐标Table 1 Specific position coordinates of sensors

图1 实验装置尺寸及布置的传感器位置示意图Fig.1 Diagram of size of experimental device and sensor location

1.2 守恒方程

MPPBuoyantEpplFoam 求解器采用了先进的算法模型，具有计算准确、稳定性高、网格适应性好等优点；同时由于其可以采用较大的CFL数进行模拟，整体计算效率较高；此外，该求解器还引进了最新开发的全速域求解算法，可以开展亚声速到超声速的工况模拟。该求解器采用多组分瞬态可压缩流动的控制方程，包括连续方程、动量方程、能量方程和组分方程。

连续方程和动量方程分别为：

式中：ρ为密度，kg/m3；t为时间，s；V为速度，m/s；p为压强，Pa；g为重力加速度，m/s2；τeff为有效剪应力，Pa。

其中，有效剪应力的表达式为：

式中：μ为混合物的分子黏度，Pa·s；μt为混合物的湍流动力黏度，Pa·s。

对于多组分流，比总焓方程形式如下：

式中：H为比总焓，J/kg；D为压力随时间变化的扩散系数，m2/s；q为热扩散通量矢量，kg/(m2·s)。

其中，比总焓的定义为：

式中：h为比焓，J/kg。

热扩散通量矢量的定义为：

式中：λ为混合物的平均导热系数，W/(m·K)；Cv为定容比热容，J/(kg·K)；Prt为湍流普朗特数；T为气体温度，K；Ns为组分的个数；hi为组分i的比焓，J/kg；ji为组分i的质量扩散通量矢量，kg/(m2·s)。

组分i的质量扩散通量矢量的定义为：

式中：Dim为组分i扩散到其他组分的平均扩散系数，计算方法参考文献[12]；Sct为湍流施密特数；Yi为组分i的质量分数。

组分i的质量分数需满足如下守恒方程：

此外，假定压强和密度满足理想气体状态方程，即：

式中：R为摩尔气体常数，J/(kmol·K)；Wi为组分i的摩尔分子质量，kg/kmol。

以上各式中，混合物的分子黏度、混合物的平均导热系数和组分i扩散到其他组分的平均扩散系数的取值与混合物的温度和组成有关，通过多组分库Mutation++计算；混合物的湍流动力黏度等湍流参数与湍流模型相关，由于考虑的实验条件受湍流的影响相对较弱，本文暂不予以考虑。

1.3 几何模型及网格划分

根据前述实验装置的尺寸和传感器的布置情况，建立几何模型并划分网格。几何模型如图2所示。

图2 几何模型Fig.2 Geometric model

根据几何模型的结构特点，本文采用结构化网格，基本单元为六面体，单元最小长度为0.5 mm，网格总数为213600个。模型泄漏口直径为0.2 cm，为提高数值模拟的准确性和计算稳定性，在泄漏口附近的X坐标方向和Y坐标方向对网格进行了加密。几何模型的计算域网格划分如图3所示。

图3 几何模型的计算域网格划分Fig.3 Mesh generation in computational domain of geometric model

1.4 模型设置

设定初始状态下模型内充满静止空气，初始压力为1 atm(1 atm=101325 Pa)(绝对压强)，初始温度为298 K。氢气从模型顶部泄漏口竖直向下注入计算域，泄漏速度为2.2 m/s，氢气注入持续时间为60 s。本模型考虑重力影响，重力方向为竖直向下，大小为9.8 m/s2。模型总计算时间为220 s，时间步长为0.0001 s，结果输出间隔为1 s。

为了探究不同位置氢气浓度随时间的变化关系并与文献[3]中的实验数据进行对比，按照表1所示的传感器具体位置坐标在模型中设置了氢气体积分数监测点。

2 实验结果分析

氢气自模型顶部泄漏后，在计算域内泄漏口处的泄漏速度切片如图4所示。

图4 泄漏口处的泄漏速度切片图Fig.4 Slice diagram of leakage velocity at leakage port

从图4可知：氢气从泄漏口竖直向下喷出，泄漏速度为2.2 m/s。喷出后扩散速度迅速减小，这主要受2个方面的影响：一方面，氢气射流受到周围静止空气的阻滞作用，发生动量交换；另一方面，氢气密度是空气密度的1/14，离开喷射口后的氢气受到较强的向上浮升力作用，从而使向下的扩散速度迅速降低。

氢气开始泄漏后，不同时间下体积分数为4%时的氢气浓度云图如图5所示，监测点处氢气体积分数随时间的变化关系如图6所示。图中：Cmax为最大氢气体积分数；th为中间的间隔时间；Δa为受限空间内上部氢气与下部氢气的浓度差，泄漏刚停止时的浓度差为12.00%，泄漏停止20 s 后的浓度差为8.03%。

从图5和图6可以看出：氢气泄漏后的整个过程可以分为注入阶段、扩散阶段和均匀混合阶段3个阶段。在注入阶段，氢气在初始动量的作用下快速在计算域中蔓延，氢气浓度迅速增加，并且在浮力作用下呈现分层状态；氢气注入结束，模型上部氢气浓度达到最大值，此时模型上下部氢气浓度差也达到最大值；此后进入扩散阶段，上下部氢气浓度差在分子扩散作用下逐渐减少；最终系统达到氢气和空气均匀混合的稳定状态，即均匀混合阶段。

图5 氢气开始泄漏后，不同时间下体积分数为4%时的氢气浓度云图Fig.5 Cloud diagram of hydrogen concentration when volume fraction is 4% at different time after hydrogen leakage

文献[3]中，传感器测得的实验装置中氢气体积分数随时间的变化关系如图7所示。图中：1、2、3 分别代表注入阶段、扩散阶段和均匀混合阶段；Δc为受限空间内上部氢气与下部氢气的浓度差。

对比图6和图7可以看出：图6的数值模拟结果与图7的实验结果曲线趋势基本一致。氢气泄漏后，由于氢气密度较低，其浓度分布在空间中表现出明显的分层效应。泄漏停止后，氢气在扩散作用下分层逐渐减弱，实验装置上部的氢气浓度水平降低，而下部的浓度水平升高，最后达到均匀混合。从定量上来看，数值模拟获得的最大氢气体积分数为12.77%，实验测得的最大氢气体积分数为13.75%，两者仅相差0.98%。泄漏停止时，数值模拟结果和实验结果均显示实验装置上部和下部的氢气浓度差为12%；20 s 后，数值模拟结果和实验结果均显示实验装置上部和下部的氢气浓度差为8%。实验结果显示，在停止泄漏后约130 s，实验装置内的氢气与空气达到均匀混合；而数值模拟结果在130 s 后，实验装置内的氢气和空气还未达到均匀混合，这是由于MPPBuoyantEpplFoam 求解器完全使用层流模型来进行扩散的计算分析，而在实际实验中很有可能受到了湍流的影响，因此氢气和空气混合的速度较快。

图6 监测点处氢气体积分数随时间的变化关系Fig.6 Variation of hydrogen volume fraction at monitoring points with time

图7 传感器测得的实验装置中氢气体积分数随时间的变化关系[3]Fig.7 Variation of hydrogen volume fraction in the experimental device with time measured by sensor[3]

通过与实验结果的对比可以看出，本文所采用的MPPBuoyantEpplFoam 求解器在模拟氢气泄漏后的流动和分层方面表现出了较高的准确性。氢气泄漏停止后，数值模拟结果中氢气的扩散速度比实验结果得到的慢，因此需要更长时间达到均匀混合。总体来说，MPPBuoyantEpplFoam 求解器的模拟准确性得到了验证，后续可以通过优化氢气扩散模型来提升模拟准确度。

3 结论

本文采用自主开发的MPPBuoyantEpplFoam求解器对受限空间内氢气泄漏扩散行为进行了数值模拟。结果显示，该求解器模拟得到的氢气浓度分布及变化情况与文献中的实验结果吻合较好，初步验证了多组分扩散功能的准确性。在本文所研究的受限空间内，氢气泄漏后表现出如下流动扩散规律：

1)氢气从模型顶部竖直向下喷入计算域后，在周围静止空气阻滞作用和浮升力作用下，扩散速度迅速减小，流动的动量控制区变短。

2)在氢气注入阶段，计算域氢气浓度分布会在浮力作用下呈现明显的分层效应，上部的氢气浓度高，下部的氢气浓度低。

3)氢气停止注入后，在分子扩散作用下，计算域上下部氢气浓度差逐渐变小，并最终达到均匀混合。

致谢：感谢哈尔滨工程大学的张慧杰博士在本文工作上给予的帮助和指导！