将天然气管道在必要地段架空敷设，有利于检修和防腐，但城市燃气管网分布在人口密集的建筑群区域，由于管道施工质量、运行维护、人为破坏及管道腐蚀等各种不可避免的因素，若未能及时检测和预防会引发天然气泄漏甚至爆炸事故^[1-3]，因此，天然气泄漏的安全性问题不容小觑。

王洪德等^[4]以某天然气储罐为例，对天然气瞬时泄漏的扩散规律进行了数值仿真，结果表明：天然气扩散距离和面积随泄漏量的增大而增大；大气越稳定，扩散的距离和范围越大；扩散距离随地面粗糙度的增大而减小。毛小虎等^[5-6]采用CFD的方法研究了天然气在二维无限大空间内泄漏扩散的危险区域，且考虑了有/无风的情况。张甫仁等^[7]采用数值模拟分析三维建筑群外空间内的天然气泄漏，将泄漏速度设为330 m/s。李朝阳等^[8]采用数值模拟分析了二维无限大空间内甲烷和硫化氢的泄漏扩散危险区域，给定了60 s内泄漏压力随时间变化的函数关系。李自力等^[9]对山地地形条件下的天然气泄漏危险区域进行数值模拟，并利用用户自定义函数（User Dedined Functiongs，UDF）修正风速在竖直方向的变化。王大庆等^[10]采用高斯模型分析天然气在有/无风条件下的危险区域。何思念等^[11]基于微分浅水方程对LNG从多种类型的储罐中泄漏并在水面扩散的过程进行了模拟。唐建峰等^[12]针对垂直喷射源连续泄漏扩散的特点，将液化天然气扩散过程分为重气扩散与被动扩散两个阶段，结合烟羽抬升计算建立了液化天然气扩散过程数学模型。朱渊等^[13]针对脱硫装置泄漏含硫天然气，大量释放致命H₂S进行了模拟和评估，定义了动态泄漏速率条件，真实反映了峰值影响时间、危害区域分布及地形影响作用。刘中良等^[14]应用热力学和气体动力学分别给出了天然气在临界状态和亚临界状态的泄漏速率计算公式，为天然气管输的安全评价提供了可靠依据。Ebrahimi等^[15]提出了一种数值方法，用来估算市区架空或埋地天然气管网的泄漏，结果表明：模拟结果的相对误差和相关值均小于5%，说明该方法具有较高的精确度。刘延雷等^[16]考虑到氢气与天然气的管道泄漏事故危险性不同，进行了数值模拟和对比，得出管道泄漏后氢气与天然气的不同泄漏扩散特性。Kurnia等^[17]应用CFD方法分析了矿井巷道内的空气流动行为和离散泄漏源，讨论了矿井巷道中多个泄漏源扩散的影响因素，提出一种控制流量和流向的新型分流器，其可将巷道内天然气累积区域的体积分数降低至安全线以下。燃气舱中气体管道泄漏一定时间后，气体体积分数在舱室上方的4个角落处聚积最多，是安装泄漏探测器的最佳位置^[18-20]。Mahmutoglu等^[21]研究了基于信号强度差的方法，用来定位水下燃气管道泄漏点的位置。

学者们在使用数值模拟方法分析天然气在建筑群之间的泄漏扩散时忽略了天然气泄漏速率随时间变化的特性，取泄漏速率为某一常数进行模拟的做法有失偏颇。在修正泄漏速率随时间变化的前提下，研究了城镇燃气管道泄漏压力和平行楼宇的几何特征（相对高度和间距）等因素对泄漏后大气流场和危险区域动态分布的影响。

1.   研究方法

1.1   流体运动控制方程

1.1.1   连续性方程

连续性方程是质量控制方程，根据流体力学的质量守恒原理建立，其一般形式为：

式中：ρ为密度，kg/m³；t为时间，s；v为流体流动的速度，m/s。

1.1.2   动量方程

动量方程又称N-S方程，是一组描述流体物质运动的方程。N-S方程的一般形式为：

对于不可压缩流体，N-S方程表示为：

式中：p为流体微团上的压力，Pa；μ为流体的动力黏度，Pa·s；f为作用在流体微团上的体积力，N。

1.2   湍流模型

采用标准k-ε两方程模型进行计算分析。标准k-ε模型采用雷诺时均法，通过求解湍动能k方程和湍流耗散率ε方程，得到k与ε的解，随后代入式（6）计算湍动黏度，进而得到雷诺应力的解。

k方程和ε方程分别表示为：

式中：k为湍动能，J；u_i为i方向的速度分量，m/s；x_i为i方向的坐标值，m；x_j为j方向的坐标值，m；μ_t为湍动黏度，Pa·s；Pr_k为湍动能普朗特数，其值为1；x_j为j方向的坐标值，m；β为热扩散系数；g_i为i方向的重力加速度，m/s²；Pr_t为湍流普朗特数；Y_M为可压湍流项中脉动扩张的贡献，J；S_k和S_ε为源项，J；ε为湍流耗散率；Pr_ε为湍流耗散率普朗特数，其值为1.3；在计算时系数C_μ=0.09，C_ε1=1.44，C_ε2=1.92，i和j的取值范围为1，2，3。

1.3   组分输运模型

天然气泄漏后与空气混合，需要利用组分输运方程计算局部每相的质量分数。由组分质量守恒定律可以得到气体扩散组分方程：

式中：ω为每一组分的质量分数；D_m为湍流扩散系数。

1.4   泄漏模型

根据开口系统能量守恒，得到式（8）：

式中：m为控制容积内的天然气质量，kg；c_V为天然气比定容热容，J/（kg·K）；T为天然气的温度，K；c_p为天然气的比定压热容，J/（kg·K）；T_e为泄漏口出口界面处气流的温度，K；c为泄漏口（出口）气流的速度，m/s；为天然气泄漏速率，kg/s。

天然气在初始泄漏阶段，由于管道压力较大，大气环境压力与管道内压力之比小于临界压力比（式9），因此天然气泄漏是处于临界状态下的泄漏；随着泄漏持续进行，管道压力逐渐降低，泄漏口处天然气的泄漏速率变为亚音速流动。因此，泄漏过程可以分为两个阶段，第一阶段为临界泄漏阶段，第二阶段为亚临界泄漏阶段。

式中：p_sur为临界压力，Pa；p为管道压力，Pa；n为天然气绝热指数，其值为1.3；β_cr为气体的临界压力比。

1.4.1   临界泄漏阶段

在临界泄漏阶段，环境压力与管道压力的关系满足式（9），泄漏口的进口处一直处于临界状态，因此，由热力学关系可知在式（8）中：

式中：R为气体常数，其值为391.78 J/（kg·℃）。

此时，由理想气体状态方程可以推导出管道中的天然气质量，由式（8）~式（10）可以得出天然气在0~t时间段内的天然气平均泄漏速率：

式中：m为天然气在0~t时间段内的平均泄漏速率，kg/s；为天然气初始泄漏速率，kg/s；m₀为管道中天然气的初始质量，kg。

1.4.2   亚临界泄漏阶段

当式（9）不能被满足，即大气环境压力与天然气管道压力之比大于临界压力比时，系统处于亚临界泄漏阶段。亚临界泄漏阶段的特点是泄漏口处气体的流动状态是亚临界的，泄漏量无法达到饱和。将式（8）结合理想气体状态方程和绝热方程，对时间求导后得到亚临界阶段泄漏速率的方程：

式中：ϕ为速度系数；A为泄漏口面积，m²；V为泄漏管道的有效容积，m³；p₀为管道初始压力，Pa；T₀为初始状态下天然气的温度，K。

1.5   模型假设与边界条件

1.5.1   模型假设

实际的天然气泄漏相当复杂，不便于理论研究，需要对其简化处理，具体假设为：①天然气是含有单一组分甲烷的理想气体，忽略泄漏过程中可能由温度引发的相变和天然气在空气中泄漏的密度变化；②泄漏后上下游阀门自动关闭，则可简化为定容泄漏问题；③泄漏源唯一，且是直径已知的一圆孔；④泄漏扩散的流态可以认为是泄漏气体与室外介质相互作用的湍流；⑤气体泄漏后只考虑其发生爆炸的可能性，不考虑其他化学反应；⑥忽略泄漏扩散过程中的热量交换与传递，即为绝热过程；⑦忽略天然气与室外介质中较大颗粒（如小水滴）之间的动量交换；⑧风向为水平方向，风向以及风速不随时间、空间而变化，楼宇外壁和地面给定速度为0的无滑移条件。

1.5.2   几何模型与边界条件

计算的流体域为30 m×35 m×50 m的长方体，其内有两栋多层建筑，A楼高12 m，B楼高18 m。GB 50352—2005《民用建筑设计通则》规定，除住宅之外的民用建筑，高度不超过24 m为单层或多层建筑，泄漏口位于两栋楼正中间，中心坐标为（-7，0，0）（图 1）。

图  1  楼宇布局物理模型示意图

下载: 全尺寸图片 幻灯片

风的进口方向与XY面垂直，风速2 m/s，XZ面为地面，其他面为出口。所研究的天然气管道的初始压力分别为0.8 MPa和1.6 MPa（代表不同工况），管径为0.66 m，天然气温度为42 ℃，平均分子量为21，环境压力为大气压，温度取25 ℃，泄漏口为直径0.2 m的圆孔，其泄漏速率随时间的推移有较大变化，由喷射速度随时间的变化曲线拟合得到质量流率与时间的关系式为m=-0.000 000 4 t³+0.000 6 t²-0.277 3 t+46.992。

多相流模型采用Mixture模型，对混合物求解动量方程、连续性方程及能量方程，对次相（甲烷）求解体积分数；湍流模型选用标准k-ε模型；采用组分传输模型模拟泄漏天然气与空气的混合；主相为空气，次相为甲烷，压力与速度的耦合方式选用SIMPLE算法，空间离散格式中，梯度离散格式选用Least Square Cell Based，压力离散格式选用PRESTO!，动量、体积分数及湍动能的离散格式选用一阶上风格式，残差收敛标准设定为0.001。

2.   结果讨论

结合现有城镇燃气管道状况以及工程应用，采用控制变量法分别研究城镇燃气管道以管道压力和平行楼宇的几何特征等为变量时天然气泄漏扩散的瞬态变化情况（表 1）。

表  1  城镇燃气管道泄漏模拟工况

下载: 导出CSV 

| 显示表格

2.1   泄漏扩散行为

一般认为，天然气爆炸极限介于0.05~0.15之间。体积分数分布图中为了标明爆炸危险范围，仅显示甲烷体积分数在0.05~0.15范围内的天然气体积分数；为展示天然气泄漏扩散形成的爆炸危险区域在空间中的具体分布，采用体积云图表示体积分数分布。由工况2条件下天然气泄漏扩散行为（图 2）可知，当管道压力为0.8 MPa时，天然气泄漏扩散的影响区域较小，两楼之间形成的旋涡较小。当天然气扩散到80 s左右时达到稳定状态，从两楼之间散开。从开始泄漏到达到稳定状态，天然气主要分布在靠近A楼背风侧的区域，因而整个扩散过程中只有A楼背风侧是危险区域。

图  2  工况2条件下城镇燃气管道泄漏CH₄体积分数分布图

下载: 全尺寸图片 幻灯片

分析工况3条件下不同时刻城镇燃气管道天然气泄漏体积分数分布情况（图 3）：天然气开始泄漏时，由于建筑物对风的阻碍，天然气从泄漏口喷射出去后四周压力平衡，因此沿泄漏口中心的竖直线呈圆锥形分布，从泄漏口向上一段距离后，天然气体积分数降低，达到爆炸危险范围，此时A、B两栋楼还未接触到天然气，均为安全区域；泄漏一段时间后（t=21.2 s），天然气沿竖直方向到达计算域顶部，B楼对局部风场的阻碍导致靠近A楼顶部的速度矢量为逆风方向，风在A、B楼之间的流场中形成一个顺时针漩涡（图 4a），泄漏天然气的主扩散中心形成一条从泄漏口至计算域顶部出口略向逆风方向偏移的扩散中心线。随着扩散的持续，天然气受风力影响，竖直向上的射流发生弯曲，向右侧扩散，从11.2 s至101.2 s，天然气呈拱形越过B楼，将B楼包围，此时整个B楼均为危险区域。

图  3  工况3条件下城镇燃气管道泄漏CH₄体积分数分布图

下载: 全尺寸图片 幻灯片

图  4  工况3在t＝21.2 s时楼宇中心的平面速度矢量图

下载: 全尺寸图片 幻灯片

随着天然气泄漏速度减缓，121.2 s后，天然气将无法越过B楼，两栋楼之间形成的回流无法散去。天然气在A楼的背风面形成副扩散中心，在A、B楼之间扩散。在121.2~221.2 s内，A、B楼之间的天然气体积分数几乎一直在爆炸极限范围内，A楼背风侧和B楼迎风侧均为危险区域。天然气扩散261.2 s后达到稳定状态，稳定后的天然气泄漏速度减缓，无法到达B楼的高度，A、B楼间流场中顺时针方向的漩涡带动更多天然气向A楼靠拢，并紧贴A楼向上扩散（图 3h、图 3i）。在XZ面上，速度矢量是由两楼间指向外部空白区域的（图 4b），泄漏口指向A楼背风方向的速度矢量导致天然气紧贴A楼分布，使得天然气从A楼顶端和两侧散开，此时A楼为危险区域。

2.2   泄漏压力对天然气扩散特性的影响

在计算过程中，分别对各工况下流体域中不同高度水平面上天然气体积分数的最大值进行监测（图 5）。

图  5  不同工况下各截面CH₄体积分数变化图

下载: 全尺寸图片 幻灯片

由工况2中天然气在不同高度截面的最大体积分数随时间的变化情况（图 5a）可知：泄漏初期，由于天然气喷射到不同高度存在时间差，不同高度的速度场分布情况不同且随时间变化，因此泄漏时间对各点天然气最大体积分数影响很大，扩散到达稳定状态后，计算域内高度10 m以下的点最终都能达到稳定状态，其天然气最大体积分数不再随时间发生变化；虽然泄漏空间内每个高度截面的天然气最大体积分数最终都能达到稳定，但因水平面高度的不同，天然气的扩散行为对该平面内天然气体积分数的影响亦不同，因此达到稳定的路径存在差别；较高位置的截面上天然气最大体积分数随时间先增大后减小，随后在波动中达到稳定；随着高度的增加，天然气被空气稀释，在各个截面上最大体积分数逐渐减小，天然气在低处体积分数大，在高处体积分数较小。

对比分析工况2和工况3（图 5b）在同一高度水平截面上的天然气最大体积分数分布，可知：①在天然气扩散过程中，天然气在某一截面的体积分数随着天然气泄漏压力的增大而增大；②天然气泄漏速率越大，其扩散达到的高度越高；③在工况2下，泄漏115 s左右，天然气在各截面上的最大体积分数达到稳定值，在工况3下，泄漏280 s左右，天然气在各截面上的最大体积分数才能达到稳定值，可见天然气泄漏速率越大，其扩散达到稳定所需的时间越长；④天然气泄漏速率越大，泄漏扩散过程中的危险区域越大。

2.3   建筑物相对高度对天然气扩散特性的影响

工况3和工况4（图 5c）的设定可以考察建筑物相对高度对天然气扩散特性的影响。在流体域中不同高度水平面上监测到天然气体积分数的最大值。工况3中A、B楼相对高度为2/3，工况4中A、B楼相对高度为1。工况4泄漏达到稳定状态后各截面天然气最大体积分数较工况3大，这意味着相对高度越大，天然气扩散达到的高度越高。B楼越高，则风速在A、B两楼之间形成的涡流速度越大，阻碍泄漏的天然气向上扩散，无法达到更高的高度。

2.4   楼距对天然气扩散特性的影响

工况1中楼距为5 m，瞬态模拟得到其各代表性时刻的爆炸危险区域（图 6）。

图  6  工况1条件下各代表性时刻CH₄体积分数分布图

下载: 全尺寸图片 幻灯片

对比工况1和工况2的危险区域，分析楼距分别为5 m和12 m时天然气体积分数分布情况及危险区域，可知：

（1）天然气泄漏扩散形成的爆炸危险区域与建筑物之间的距离有密切关系。楼距越大，天然气扩散速度越快，天然气越容易散开；反之越容易积聚在两楼之间。

（2）楼距越大，天然气受到两楼间顺时针旋涡的影响越大，天然气扩散达到的高度越高，在扩散过程中，天然气分布较广，危险区域较大。

（3）楼距为5 m时，天然气扩散90 s左右达到稳定状态，但此时A、B两楼之间依然是危险区域，存在较大的安全隐患；楼距为12 m时，天然气扩散120 s左右达到稳定状态，稳定后的天然气残余量明显较小，长期来看，天然气扩散稳定后的残余量越少，说明天然气扩散越快，安全性相对越好。因此适当增加楼距有利于加快天然气扩散，减少天然气泄漏后的安全隐患。故GB 50028—2006《城镇燃气设计规范》中规定次高压燃气管道距建筑物外墙面不应小于3 m。

3.   结论

采用瞬态数值模拟方法，考虑天然气泄漏速率随时间的动态变化，针对平行楼宇间天然气管道泄漏扩散形成的危险区域进行分析，得到以下结论：

（1）天然气扩散过程中遇到建筑物，会在其背风向形成副扩散中心向四周扩散；随着距离地面高度的增加，天然气最大体积分数逐渐降低，且随着泄漏持续达到某一稳定值；楼距较大的工况下A楼背风侧天然气均呈现贴壁流动，因此楼宇内的居民处于危险区域，B楼的危险是暂时的，而楼距小的工况下，A、B两楼均为危险区域，若天然气扩散进入房间，有可能引发爆炸事故。

（2）随着管道压力的增大，天然气泄漏扩散的影响区域和危险区域均增大，因此可以通过适当减小管输压力从而缩小危险区域；相邻楼宇的相对高度对天然气危险区域的分布有较大影响，有风的情况下，若泄漏口位于较高楼宇的上风向，则上游为危险区域，下游较安全；增大楼距虽然使得天然气扩散范围更广，扩散过程中的危险区域更大，但可以将天然气更快地散布出去，安全性更好。