基于FLACS的天然气泄漏扩散数值模拟研究*

罗振敏 胡腾 李俊 王涛 高有军 杨勇 王登飞

（1.西安科技大学安全科学与工程学院，西安 710054；2.陕西省工业过程安全与应急救援工程技术研究中心，西安 710054；3.陕西省煤火灾害防治重点实验室，西安 710054；4.榆林榆川天然气有限责任公司，陕西 榆林 719099）

0 引言

近年来我国天然气泄漏着火爆炸事故屡有发生。2021年10月18日，河北省邯郸市发生天然气泄漏事故，造成3人窒息死亡。2021年10月21日，辽宁省沈阳市和平区一家饭店发生天然气泄漏爆炸事故，造成3人死亡、30多人受伤。2021年6月13日，湖北省十堰市发生重大天然气爆炸事故，致26人死亡。2021年10月24日，辽宁省瓦房店市一居民楼发生燃气闪爆事故，造成2人死亡、7人受伤。因此，研究天然气在各种工况下的泄漏扩散规律及浓度、速度分布规律对人员提前采取防护措施和制定火灾爆炸预防方案具有重要意义。数值模拟研究具有可操作性强、成本低、应用面广等优势［1］。

国内外大量学者对天然气泄漏扩散进行了数值模拟研究，其中基于CFD模型进行研究的有：王春园［2］在FLUENT中对管道正常情况和泄漏情况分别进行了仿真，得到了燃气管道泄漏后管道首末两端气体压力和流量的变化规律；YUE C J等［3］利用计算流体动力学软件FLACS，对徐州某LNG储罐泄漏可能造成的多种灾害进行了模拟分析；王文和等［4］以重庆某无人值守井站为实际场景，采用ANSYS软件对含硫天然气泄漏扩散过程进行模拟。基于ALOHA软件进行研究的有：彭琳等［5］为提高地面天然气管道泄漏扩散范围预测的精度，基于ALOHA事故后果模拟分析和多元回归预测方法建立地面天然气泄漏扩散范围预测模型。基于高斯烟羽模型、内嵌UDM模型的DNV（挪威船级社）PHAST、SAFETY软件进行研究的有：何宁辉等［6］通过DNV公司的SAFETI软件对于某储配厂的天然气罐区进行了泄漏后果的模拟和定量风险分析。基于开源CFD软件Open-FOAM进行研究的有：周理［7］利用开源CFD计算程序包OpenFOAM对激波管内的高压燃气泄漏自燃现象进行了数值模拟。

以上研究都未对气云速度矢量在不同方向的分量分布规律进行分析。本文基于FLACS软件，首先对文献［7-8］中的物理实验得出的天然气泄漏扩散浓度与速度结果进行数值模拟可靠性验证，然后对某天然气储罐区进行泄漏扩散模拟，得出天然气在一定工况下泄漏扩散的浓度、速度矢量分量、等效可燃气云体积Q9变化规律。以期对处置天然气储配气站场泄漏扩散事故具有指导意义。

1 模型建立与场景设置

本文模拟对象为某天然气储配站，如图1所示，罐区共有8个球形钢制储罐（编号L1—L8）和3间控制室和厂房，每个球罐体积约为2 000 m3（球罐直径约16 m，底部距地面5.0 m），球罐与球罐之间布置有扶梯和管道。

图1 天然气储罐区CFD几何模型

在FLACS软件CASD建立相应的天然气储罐区CFD几何模型。假设正北方向为+Y，计算区域坐标为：X（-5～248 m）；Y（-20～115m）；Z（0～35 m）；图左下角为坐标原点。对CFD模型划分网格，并对泄漏源附近网格细化，在主要研究区域之外，对网格拉伸，在细网格和粗网格之间进行平滑操作，最终划分网格982 688个。然后对网格模型进行孔隙度计算。

天然气在泄漏扩散过程中遵循质量守恒方程、能量守恒方程、动量守恒方程。

2 浓度与速度结果可靠性验证

本文分别对文献［7-8］中的小尺寸物理实验进行浓度及速度结果验证。在FLACS中建立与文献描述完全相同的CFD几何模型，按文献所述进行泄漏源位置、泄漏方向、泄漏孔径及其他初始条件的设置，对监测点进行适当删减，只保留文献［7］模型中心浓度监测点和文献［8］4号速度测点。浓度结果及速度结果的实验值与模拟值对比如图2—图3所示：

图2 浓度结果对比

图3 速度结果对比

由图2—图3可得利用FLACS的模拟结果与实验数据的变化趋势相近，但局部误差较大，分析误差原因为：在现场实验过程中大气条件发生一定变化，而数值模拟过程中假设大气条件不改变。此外，文献中部分初始条件没有详细列出，可能与模拟中的初始条件有差异，最终导致产生误差。图2实验与模拟数据平均相对误差为15.36%，图3实验与模拟数据平均相对误差为11.45%，误差均在20%以内，符合工程允许误差范围之内。综上所述，FLACS数值模拟得出的浓度与速度结果相对可靠。

3 数值模拟结果分析

本文模拟L2储罐底部发生泄漏，泄漏源坐标为（45 m，50 m，5 m），泄漏孔径40 mm，泄漏方向为+Y（正北）。环境温度20℃，大气压101 325Pa，地面粗糙度0.01，大气稳定度等级F（由环境风速、云覆盖面积等天气情况决定），泄漏时间为60 s，泄漏停止后持续扩散时间为30 s，具体工况见表1。

表1 天然气泄漏扩散工况

3.1 风速对气云扩散的影响

图4—图7分别为case1、case2、case3、case4条件下泄漏后的天然气扩散达到稳定状态时Z=5 m的XY切面（泄漏源所在平面）天然气体积分数（1.5%～10%）分布图、天然气速度矢量（VVEC）在V方向上的分量（-35～35 m/s）分布图、天然气速度矢量（VVEC）在U方向上的分量（-10～10 m/s）分布图、等效可燃气云体积Q9随时间变化图。

图7 不同风速Q9随时间变化值

由图4可得风速越大，天然气云扩散到达的最远距离越近，分别为55、45、40、37m。无风状态时，浓度沿泄漏源中轴线左右对称分布，存在横向风力时，风速越大，射流与泄漏源中轴线倾斜角度越大。射流核心区域浓度远高于边界层，这是因为天然气喷射出后，其气云边界层不断与大气进行质量、动量、能量交换，导致边界层不断卷吸空气，浓度变低，风速越大，边界层卷吸空气程度越强。同时，由于风力的横向输送作用，风速越大，气云团横向输送距离越远，射流中轴线倾斜角度越大。总的来看，4种风速条件下射流中轴线倾斜角度均较小，这是因为模拟泄漏速率较大，动量主导气云团运动轨迹［9］。

图4 不同风速XY切面气云浓度

由图5可得，速度矢量（VVEC）在V方向的分量分布图与浓度分布图类似，其边界层值较小，射流中心区域值较大，这是因为气云团边界层不断与大气进行动量交换，导致速度衰减。又因为风力的横向输送作用，速度云图与泄漏源中轴线形成一定角度，风速越大，角度越大。此外，同种风速情况下泄漏源远端速度云图倾斜角度大于近端，这是因为天然气刚离开泄漏孔时速度远大于风速，风速对其喷射影响很小，随着径向距离增大，天然气速度逐渐衰减，大气湍流的作用增强，天然气速度等值线向右倾斜的趋势越来越明显［10］。

图5 不同风速XY切面气云速度矢量在V方向的分量分布

由图6可得，速度矢量（VVEC）在U方向的分量分布图与浓度云图不同，无风时呈现“蝴蝶”形状左右对称分布，这是因为射流核心区域卷吸边界层，导致对称轴右侧的天然气产生向-X方向的速度，对称轴左侧的天然气产生+X方向的速度。风速逐渐增大时，监测区域平面的气云速度矢量在U方向的分量分布逐渐不规则，反映出湍流程度增强，但总体在射流中轴线左侧的气云速度大于右侧气云速度，这也是由于射流核心区域卷吸边界层，导致射流轴线右侧的天然气产生向-X方向的速度，对称轴左侧的天然气产生+X方向的速度。

图6 不同风速XY切面气云速度矢量在U方向的分量

FLACS中用等效可燃气云体积Q9来表示阻塞率低，通风良好的开放场景与真实气云等效化学计量比的均匀气云体积，m3。Q9的计算公式为［9］：

式中，V为可燃气体体积；BV为层流燃烧速度；E为在空气中恒压燃烧所产生的体积膨胀。

由图7可得不同风速条件下等效可燃气云体积Q9随时间变化趋势基本一致。随着泄漏开始，Q9在极短时间内从0迅速增大，随后增速放缓，最终达到最大值，此时泄漏仍在持续，但Q9只发生小幅度波动，基本保持最大值不变。60 s时，随着泄漏停止，Q9值迅速下降，在2 s之内从最大值降为0。从纵向相比可以看出，在相同时刻，风速越大，监测区域内的等效可燃气云体积越小，这是因为随着风速的增大，监测区域内气体湍流程度增强，天然气卷吸空气程度增强，其浓度稀释程度增强，导致Q9变小。此外，4种风速条件下Q9值趋于稳定后，风速越大，Q9值波动越大，这也是由于风速越大，气体湍流程度越强的原因造成的。

3.2 泄漏速率对气云扩散的影响

图8—图11分别为case1、case5、case6条件下泄漏后的天然气扩散达到稳定状态时Z=5 m的XY切面（泄漏源所在平面）天然气体积分数（1.5%～10%）分布图、天然气速度矢量（VVEC）在V方向上的分量（-35～35 m/s）分布图、天然气速度矢量在U方向上的分量（-10～10 m/s）分布图、等效可燃气云体积Q9随时间变化图。

由图8可得，气云图沿泄漏源中轴线左右对称分布，随着泄漏速率增大，天然气云团扩散到达的最远距离逐渐增大，分别为37、55、95 m。这是因为当泄漏孔径一定（40 mm）时，泄漏速率越大，天然气刚喷射出时初始动能越大，在与大气进行质量、动量、能量交换过程中，其动能衰减到0所需时间越长，导致最终扩散的最远距离越远。而且，泄漏速率越大，相同时间内所喷射出的气云体积越大，导致气云覆盖范围越大。

图8 不同泄漏速率XY切面气云浓度

由图9可得，与浓度分布云图类似，速度矢量（VVEC）在V方向的分量分布（-35～35 m/s）云图沿泄漏源中轴线左右对称分布。随着泄漏速率增大，V方向上的速度分量到达最远距离逐渐变大，分别为68、89、101 m，这是因为泄漏孔径一定（40 mm）时，泄漏速率越大，刚喷射出的天然气动量越大，在气云卷吸大气的过程中，其速度衰减到0的过程越慢，导致速度边界层到达距离越远。

图9 不同泄漏速率XY切面气云速度矢量在V方向的分量

由图10得，速度矢量（VVEC）在U方向的分量（-10～10 m/s）分布图均呈“蝴蝶”形状左右对称分布，这是因为泄漏孔处形成负压区，射流不断与大气掺混造成的。随着泄漏速率的增大，速度云图分布面积增大。这是因为泄漏速率越大，射流核心区域对大气的卷吸作用越强造成的。除此之外，计算域边界出现面积较大的速度云图，这是由于射流喷射到最远端时，沿核心区域对称轴形成“涡团”，继续卷吸大气造成的，泄漏速率越大，“涡团”尺寸越大［11］。

图10 不同泄漏速率XY切面气云速度矢量在U方向的分量分布

由图11可得3种泄漏速率Q9变化趋势基本相同，但从纵向相比可以看出，在相同时刻，泄漏速率越大，监测区域内的等效可燃气云体积越大，这是因为泄漏速率越大，单位时间内扩散到浓度监测区域的天然气体积越大，导致Q9值越大。

图11 不同泄漏速率Q9随时间变化值

4 结论

1）通过对天然气泄漏扩散浓度分布和速度分布实验结果的数值模拟验证，证明FLACS数值模拟在研究天然气泄漏扩散浓度及速度分布规律具有一定可靠性；

2）横向风速越大，天然气扩散最远距离越近，浓度云图向下风向倾斜的趋势越明显；速度矢量在V方向的分量分布与浓度分布类似，风速增大会使速度等值线向右倾斜的趋势明显；随着横向风速增大，速度矢量在U方向的分量分布越不规则，但总体呈沿射流中轴线左侧的气云速度大于右侧气云速度；等效可燃气云体积Q9趋于稳定后值随横向风速增大而减小，且波动幅度增大；

3）泄漏速率越大，天然气刚喷射出时初始动能越大，速度衰减到0所需时间越长，扩散距离越远，稳定后形成气云覆盖范围越大，速度矢量在V方向的分量边界层到达距离也越远；由于射流核心区域的卷吸大气作用，射流中轴线左侧形成天然气产生+X方向速度，右侧产生-X方向速度，速度矢量在U方向的分量分布云图在泄漏源附近呈“蝴蝶”状分布，泄漏速率越大，分布面积越大；射流末端出现“涡团”，泄漏速率越大，“涡团”尺寸越大；在相同时刻，泄漏速率越大，监测区域内的等效可燃气云体积Q9越大。