燃气管网是城市的关键基础设施之一，现代城市的正常运行对燃气供应有极强的依赖性。但是，燃气管网泄漏严重威胁公共安全，一旦发生爆炸极易造成人员伤亡和财产损失。2014年3月3日中国东莞发生燃气爆炸事故，导致1死31伤。2014年8月1日，中国台湾高雄市多条街道因可燃气体泄漏发生连环爆炸，造成32人遇难，321人受伤。因此，城市燃气管网的安全运行对保障人民生命财产安全和社会正常运转具有重大的意义。城市燃气管网大多是埋地敷设，泄漏的燃气需要先在地下土壤环境中扩散，然后才通过地面进入街道在大气环境中扩散。研究燃气在土壤环境和大气环境中的耦合扩散规律对燃气泄漏的监测预警和应急处置有重要意义。

对于燃气在土壤中的扩散规律，国内外学者从实验、数值模拟和理论分析3个方面开展了相关研究。孙立国等^[1]基于德国水和燃气协会 (DVGW) 的系列实验对天然气在地下扩散的影响因素进行了分析，包括温度、湿度、路面是否覆盖沥青或水泥和路面是否有裂隙或检查井等，并指出正确掌握燃气的泄漏量、泄漏速度、扩散范围和浓度分布等关键信息对现场救援非常重要。Praagman等^[2]进行了地下天然气运移的多相流数值模拟和小尺寸实验研究，指出土壤的均匀性、饱和度、温度和管线压力等都是研究天然气扩散范围时需要考虑的因素。Okamoto等^[3]针对埋地天然气管道泄漏扩散进行了全尺寸实验研究，得到地下甲烷浓度分布随时间的变化情况，并通过实验与三维数值模拟的对比，验证了Darcy定律与Fick扩散定律在全尺寸地下管道泄漏扩散问题中的适用性。谢昱姝等^[4]采用自制的全尺度气体泄漏实验系统，对泄漏后天然气在土壤中的扩散过程进行研究，指出埋地管道泄漏后天然气在土壤中的扩散过程可以分为孕育、陡然增长、缓慢增长和稳定4个阶段。熊兆洪等^[5]基于单相渗流和高度非Darcy渗流理论建立了埋地天然气管道微小泄漏的数值模型，使用有限体积法求解控制方程，并和小尺度沙槽空气渗流实验的结果进行了对比。晏玉婷等^[6]使用FLUENT软件对埋地中压天然气管道的泄漏进行了数值模拟，并通过全尺寸实验比较了低压天然气管道在4种泄漏方向下的浓度分布情况，指出泄漏方向对管道上方的甲烷浓度有显著的影响^[7]。张鹏等^[8]针对埋地天然气管道的微小泄漏，通过忽略土壤中甲烷对流扩散方程的对流项得到了甲烷在土壤中扩散的解析表达式，并通过参数的敏感性分析指出泄漏孔径、土壤含水率和土壤种类等均是影响扩散过程的重要因素。

在土壤与大气耦合的燃气扩散方面，Esposito等^[9]通过小尺寸的风洞实验研究了大气边界层对土壤中甲烷扩散规律的影响。Botros等^[10]分别为土壤和大气中的燃气扩散构建了解析解，并通过地面燃气扩散通量进行了土壤大气扩散的耦合。Parvini等^[11]对土壤中的燃气扩散进行了数值模拟，并将计算出的地面扩散通量作为地面上扩散、爆炸等解析模型的输入条件，对埋地管道泄漏进行了风险评估。程猛猛等^[12]使用FLUENT软件自带的基于阻力源项法的多孔介质模型对不同的天然气泄漏情景进行了土壤、大气耦合的数值模拟。

目前，在燃气泄漏数值模拟中考虑土壤、大气耦合的研究较少，且多采用附加阻力源项的方法来表征土壤，很少有在土壤和大气2个区域采用不同控制方程的。本文面向城市地下燃气泄漏探测和评估的需求，针对土壤和大气中燃气扩散的物理规律分别构建控制方程，并基于开源计算流体力学工具OpenFOAM进行二次开发，对控制方程进行数值求解，开展土壤、大气耦合的燃气扩散数值模拟研究，重点分析了地面对流传质系数对冒出地面的甲烷通量和街道内甲烷浓度测量值的影响。

1 物理模型 1.1 土壤扩散物理模型

考虑与文[3]中实验类似的单相微小泄漏情形，此时土壤中甲烷的扩散规律可以用Darcy定律和Fick定律来描述。对于符合Darcy定律的土壤渗流，其气体渗流速度v_g与土壤的渗透率k，流体的动力黏度μ，压力梯度$\nabla p$，密度ρ和重力加速度g有关，表达式为

$ {\mathit{\boldsymbol{v}}_{\rm{g}}} = \frac{k}{\mu }\left( { - \nabla p + \rho \mathit{\boldsymbol{g}}} \right). $

(1)

对于混合气体，动力黏度μ依据文[13]获得

$ \mu = \sum\limits_i {\frac{{{\mu _i}}}{{1 + \sum\limits_{j \ne i} {{\phi _{ij}}{x_j}/{x_i}} }}} , $

(2)

$ {\phi _{ij}} = \frac{{{{\left\{ {1 + {{({\mu _i}/{\mu _j})}^{\frac{1}{2}}}{{({M_j}/{M_i})}^{\frac{1}{4}}}} \right\}}^2}}}{{2\sqrt {2(1 + {M_i}/{M_j})} }}. $

(3)

其中：μ_i是第i个气体组分的动力黏度；x_i是第i个气体组分的摩尔分数；M_i是第i个气体组分的摩尔质量。

用组分质量分数C_i作为变量的组分输运方程如下：

$ \frac{\partial }{{\partial t}}(\varepsilon \rho {C_i}) + \nabla \cdot(\rho {C_i}{\mathit{\boldsymbol{v}}_{\rm{g}}}) = \nabla \cdot(\rho D\;\nabla {C_i}) + {S_i}. $

(4)

其中：ε是土壤的孔隙度；D是气体在土壤中的有效扩散系数；S_i是第i个气体组分的质量源项。

连续性方程如下：

$ \frac{\partial }{{\partial t}}\left( {\varepsilon \rho } \right) + \nabla \cdot(\rho {\mathit{\boldsymbol{v}}_{\rm{g}}}) = \sum\limits_i {{S_i}} . $

(5)

理想气体状态方程如下：

$ p = \sum\limits_i {\frac{{{\rho _i}RT}}{{{M_i}}}} = \sum\limits_i {\frac{{{C_i}\rho RT}}{{{M_i}}}} = \frac{{\rho RT}}{M}. $

(6)

其中：ρ_i是第i个气体组分的密度；R是气体常数，取值为8.31；T (K) 是温度。

将式 (1) 代入式 (5)，并删掉非稳态项，得到关于压力的方程，即

$ \nabla \cdot\left( {\rho \frac{k}{\mu }(\nabla p)} \right) = - \sum\limits_i {{S_i}} + \nabla \cdot\left( {\rho \frac{k}{\mu }\rho \mathit{\boldsymbol{g}}} \right). $

(7)

1.2 大气扩散物理模型

采用标准k-ε模型^[14]，求解大气中的气体扩散过程，控制方程如下 (注意，这里变量的符号遵循CFD领域的惯例，与节1.1中的符号是独立的)：

$ \frac{{\partial \overline {{u_i}} }}{{\partial {x_i}}} = 0, $

(8)

$ \frac{{{\rm{D}}\overline {{u_i}} }}{{{\rm{D}}t}} = - \frac{1}{\rho }\frac{{\partial \bar p}}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_j}}}\left( {\nu \frac{{\partial \overline {{u_i}} }}{{\partial {x_j}}} - \overline {{{u'}_i}{{u'}_j}} } \right), $

(9)

$ \frac{{{\rm{D}}k}}{{{\rm{D}}t}} = {P_k} + \frac{\partial }{{\partial {x_j}}}\left( {\frac{{{\nu _t}}}{{{\sigma _k}}}\frac{{\partial k}}{{\partial {x_j}}}} \right) - \varepsilon , $

(10)

$ \frac{{{\rm{D}}\varepsilon }}{{{\rm{D}}t}} = \frac{\varepsilon }{k}({C_{\varepsilon 1}}{P_k} - {C_{\varepsilon 2}}\varepsilon ) + \frac{\partial }{{\partial {x_j}}}\left( {\frac{{{\nu _t}}}{{{\sigma _\varepsilon }}}\frac{{\partial \varepsilon }}{{\partial {x_j}}}} \right), $

(11)

$ \overline {{{u'}_i}{{u'}_j}} = {\nu _t}\left( {\frac{{\partial \overline {{u_i}} }}{{\partial {x_j}}} + \frac{{\partial \overline {{u_j}} }}{{\partial {x_i}}}} \right) - \frac{2}{3}{\delta _{ij}}k, $

(12)

$ {\nu _t} = {C_\mu }\frac{{{k^2}}}{\varepsilon }, $

(13)

$ {P_k} = - \overline {{{u'}_i}{{u'}_j}} \frac{{\partial \overline {{u_i}} }}{{\partial {x_j}}}, $

(14)

$ \frac{{{\rm{D}}C}}{{{\rm{D}}t}} = - \overline {{u_i}} \frac{{\partial C}}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_i}}}\left[ {D + \left( {\frac{{{\nu _t}}}{{S{c_t}}}\frac{{\partial C}}{{\partial {x_i}}}} \right)} \right] + {S_C}. $

(15)

其中：${\overline {{u_i}} }$是时均速度；u′_i是脉动速度；ρ=1.205是空气的密度；p是时均压力；v=1.5×10^-5是空气的运动黏度；ν_t是涡黏性系数；k是湍动能；ε是湍流耗散率；C_μ=0.09, C_ε1=1.44, C_ε2=1.92, σ_k=1.0和σ_ε=1.3是湍流模型的常数；C是组分质量浓度；D=1.5×10^-5是扩散系数；Sc_t=0.7是湍流Schmidt数。

2 计算工况与计算方法 2.1 土壤计算区域和边界条件

参照文[3]，选取图 1所示的二维区域和边界条件来计算甲烷在土壤中的扩散。计算域长30 m，深2.9 m，其中黑色圆圈为甲烷泄漏点。土壤温度取25 ℃。甲烷动力黏度取11.067 μPa· s，空气动力黏度取18.448 μPa·s。甲烷泄漏速率为14.8 mg/s。土壤孔隙率为0.17，渗透系数为2.0×10^-11 m²，有效扩散系数为0.66×10^-6 m²/s。地面边界条件中系数H取5.0/m，侧面假设为充分发展边界，底面为地下水，即为不透气边界。

注意这里的系数H并不是对流传质系数，依据传质理论，在地面处甲烷传质边界条件的一种更为完整的写法是

$ \begin{array}{c} \rho {v_n}{C_{{\rm{2,soil}}}} - \rho D\frac{{\partial {C_2}}}{{\partial n}} = \rho {h_{\rm{m}}}({C_{{\rm{2,air}}}} - {C_{{\rm{2,ref}}}}),\\ {C_{{\rm{2,soil}}}} = k{C_{{\rm{2,air}}}}. \end{array} $

其中：v_n是地面处的法向速度；C_{2, air}是地面空气侧的甲烷质量分数；C_{2, soil}是地面土壤侧的甲烷质量分数；C_{2, ref}是空气中的参考甲烷质量分数，在CFD计算中一般为第一层网格处的组分质量分数；h_m是空气和地面之间的对流传质系数；k是土壤和空气之间甲烷的分配系数。系数H可以用上述变量进行表达：

$ H = \frac{{\frac{{ - \partial {C_2}}}{{\partial n}}}}{{{C_{{\rm{2,soil}}}}}} = \frac{{{h_{\rm{m}}}(1/k - {C_{{\rm{2,ref}}}}/{C_{{\rm{2,soil}}}}) - {v_{\rm{n}}}}}{D}. $

由此可见，系数H只是规定了甲烷质量分数的梯度与甲烷质量分数的比值，它的真实取值与很多因素有关，但是增大H在一定程度上相当于增大了对流传质系数h_m，在后续的计算中以此为依据研究对流传质系数对地下扩散过程的影响。

2.2 大气计算区域和边界条件

为方便研究泄漏后的燃气在街道中的扩散规律，选取一个风洞实验的几何模型作为大气计算区域^[15-16]。该二维计算区域包含五个街道峡谷见图 2，比例尺为1:500^[15]，高宽比均为1。燃气泄漏发生在字母S标记的街道里，具体位置和尺寸在后面说明。入口处按照风洞实验的结果给定速度、湍动能和湍流耗散率的廓线，上边界取为对称边界，出口取为充分发展边界。

2.3 计算方法

使用开源CFD工具包OpenFOAM和SIMPLE算法^[17]求解式 (1)-(15)。OpenFOAM是采用C++语言编写的开源连续介质力学数值计算库，包含网格划分、构建求解器和后处理等功能。用户可以针对研究的具体问题，利用OpenFOAM和其开源社区提供的“基础设施”自定义求解器，在使用商业软件和自己编制程序之间提供了灵活的选择余地。因此，本文选择OpenFOAM作为数值计算的工具。在OpenFOAM中构建了多组分Darcy流动求解器和组分输运求解器，分别求解甲烷在土壤和大气中的扩散，两部分通过地面的甲烷通量进行耦合。所有变量的收敛准则均为残差小于10^-6。

3 结果与讨论 3.1 土壤扩散计算结果 3.1.1 无泄漏情形

当没有泄漏发生时，参考式 (1)，由于重力作用，压强随高度的增加而减小，呈现出规则的层状分布，见图 3。

3.1.2 有泄漏情形

当在 (0, 0) 发生14.8 mg/s的甲烷泄漏时，达到稳态后的压强分布如图 4，甲烷质量分数分布如图 5。可见泄漏口附近压强增大。甲烷的影响半径大概为2 m，与文[3]中的结果相近。

3.1.3 增大泄漏量

将泄漏量增大一倍，稳态时的压强分布和甲烷质量分数分布见图 6和7。可见压强分布更加偏离层状分布，甲烷的影响范围有所扩大，中心高浓度区域也明显增大。

3.1.4 增大对流传质系数

将系数H由5.0/m增大到50.0/m，得到的甲烷质量分数分布如图 8。总体看来与图 5相比没有明显的差异，但是在地面附近甲烷质量分数降低了，这是因为系数H是甲烷质量分数的梯度与甲烷质量分数的比值。系数H的影响在节3.2中进一步讨论。

3.2 地面传质通量

图 9展示了系数H对地面甲烷通量的影响，甲烷的总传质通量、扩散传质通量和对流传质通量均呈现类似钟形的分布，总的传质通量几乎不受系数H的影响。H影响的是对流传质通量和扩散传质通量的比例。当H增大时，扩散作用对地面甲烷通量的贡献增加，对流作用对地面甲烷通量的作用减小。

3.3 大气扩散计算结果

将土壤扩散计算得到的地面甲烷总传质通量作为源项，计算甲烷在大气中的扩散。如图 10所示，将地面甲烷通量分布在上风处临近建筑物的5 m宽的路面上，这与天然气管道在实际中的埋设位置是相符的。将实际通量分布和均匀的通量分布进行对比，并设置4条质量浓度观测线，以分析不同位置监测到的甲烷质量浓度对地面甲烷通量分布方式的敏感性。图中数字表示缩比模型中相对于地面的高度 (m)。

图 11对均匀通量分布和实际通量分布的计算结果进行了对比，甲烷质量浓度的数值较大是因为计算是在缩比模型中进行的。4条观测线均按照逆时针方向对甲烷质量浓度进行采样。在图 11a中，距离大代表接近地面，可见仅在地面附近两种通量分布造成的甲烷质量浓度有较大区别。在图 11b中，距离小表示街道的左侧，可见在接近泄漏源的街道左侧得到的甲烷质量浓度对通量分布比较敏感。在图 11c中，距离小代表接近地面，由于观测线C整体都在街道右侧，其质量浓度分布对通量分布不敏感。图 11d和11b反映的现象类似，由于观测线D在观测线B的上方，故其质量浓度分布对地面通量分布的敏感性有所降低。而观测线D的高度对应实际中的1.5 m，正是车载式甲烷探测器的一个有代表性的安装高度。

3.4 讨论

由上述计算结果可以知道，在一些条件下地面通量的具体分布对车载探测的影响并不明显，尤其是行车路径不在泄漏点正上方时。这使得通过车载测量值对城市燃气泄漏进行快速且较为准确的评估成为可能。事实上，出于安全和环保的考虑，基于车载式气体探测技术的天然气泄漏量估计已经成为美国等发达国家的一个研究热点，主要关注点包括采样策略与估计算法的设计、行车轨迹的规划和对估计结果不确定性的评估等^[18-20]。结合国外燃气泄漏监测、估计实践和本文计算结果，在有条件的情况下，以不同行车方向多次通过相同路段有助于提高浓度监测的可靠性。另外，如能搭载风速仪等气象监测设备将有助于通过甲烷浓度监测值反推燃气的泄漏方位和泄漏速率。

本文在构建土壤中的甲烷对流扩散模型时假设流动是等温、单相的Darcy流动，且甲烷的扩散符合Fick定律，在本文研究的微小泄漏情景中这些假设是合理的^[3]。在构建大气中的甲烷对流扩散模型时假设流动是等温的，且没有考虑甲烷密度小于空气而产生的浮力效应，这在风速较高时是合理的^[21]，在风速非常低时有必要进一步考虑甲烷的浮力效应。另外本文没有考虑复杂土壤和路面环境等非理想的情景，可在未来选取典型情景开展进一步研究。

4 结论

本文面向城市地下燃气管网泄漏探测和评估的需求，针对土壤和大气中燃气扩散规律的不同，在OpenFOAM中对不同的控制方程分别进行了求解，并通过地面的传质通量将土壤和大气中的扩散过程耦合起来。

通过数值模拟研究发现，地面的对流传质系数主要影响地面对流传质通量和扩散传质通量的比例，而对总的传质通量影响很小。只有监测点在泄漏源正上方时，测得的质量浓度值才对地面通量的具体分布比较敏感，但是此时地面传质通量分布造成的质量浓度分布差别依然很小。这些结果表明，相比土壤环境中比较稳定的燃气扩散过程，大气环境中瞬息万变的燃气扩散过程主导了实际应用中质量浓度监测的结果。下一步将通过实验研究对本文结论作进一步的检验和验证。