1. 引言

气体扩散模型[1]通常可分为两大类：一类基于数学计算，另一类基于历史数据建模。数学计算类模型的典型代表包括高斯扩散模型[2]和计算流体力学(CFD)模型。高斯扩散模型以简洁的数学公式进行计算，具有操作简便、时间成本低等优势，但其适用范围主要局限于平坦、无障碍的环境，对于复杂地形的气体扩散预测则存在精度不足的问题[3]。相比之下，CFD模型能够灵活应对复杂的几何形状，在模拟城市建筑等复杂地形中的气体扩散时表现出较高的精度，因此被广泛应用于城市环境中有毒气体扩散预测[4]。沈锴欣等[5]通过CFD模拟分析了化工多米诺事故中火灾与有毒气体泄漏的耦合效应，研究表明这种耦合作用显著加剧了事故后果的严重性。王曌文等[6]利用CFD模拟研究了双燃料动力船舶透气桅释放气体的扩散特性，分析了透气桅设计高度对气体扩散的影响，为船舶安全设计提供了科学依据。王一昊等[7]基于CFD仿真对LNG加气站的泄漏场景进行了模拟，评估了现有可燃气体探测器的覆盖率，并提出优化方案以提升覆盖率水平，从而为泄漏早期预警提供了技术支持。Tian等[8]建立了高压氢气泄漏与扩散过程的CFD模型，该模型综合考虑了泄漏过程中的热交换与真实气体状态方程，能够有效预测高压氢气泄漏行为及其对周围环境的影响。

在进行实际工程问题的数值模拟时，首先需要确保数值解的收敛性，即数值解是存在且稳定的，其次应确保数值解与网格数量无关。Freitas [9]对数值模拟的不确定性评价方法进行了系统综述，指出数值计算的收敛解必须与网格数量无关，以保证网格独立性或收敛性。Manna [10]等研究表明，对于稳态边界值问题，CFD数值误差主要来源于迭代收敛误差和网格收敛误差，且当两个连续网格的GCI值小于3%时，可视为网格已达到收敛要求。

本研究以南宁市新邕路综合管廊燃气舱防火分区9为研究对象，参考了管廊建造图纸建立三维模型，基于ANSYS Fluent软件进行了地下综合管廊燃气舱内管道泄漏的流体力学仿真模拟。研究重点包括模型验证和网格独立性分析，为后续其他气体扩散仿真模型的构建与研究提供了理论支持与方法指导。

2. 模型与计算

2.1. 物理模型

本研究针对南宁市新邕路综合管廊中燃气舱防火分区9内的天然气管道泄漏工况进行了建模分析。依据管廊建设图纸及天然气管道的实际铺设情况，构建了一个燃气舱物理模型，模型尺寸为长度180 m、宽度1.9 m、高度3.5 m。燃气舱内天然气管道直径为0.32 m，进风口和出风口均为边长0.74 m的矩形截面。泄漏模拟选取典型工况，假设泄漏口直径为0.2 m，具体模型尺寸及结构特征如图1所示。

Figure 1. Schematic diagram of pipe gallery model

图1. 管廊模型示意图

2.2. 数学模型

a) 连续性方程

$\text{div}\left(\rho u\right)=0$ (1)

式中： $\rho$ 为密度，kg/m³； $u$ 为流体速度矢量，m/s。

b) 动量方程

$\frac{\partial }{\partial t}\left(\rho u\right)+\nabla \left(\rho uu\right)=-\nabla p+\nabla \left(\tau \right)+\rho \text{g}+F$ (2)

式中： $p$ 为压力，Pa； $F$ 为外力矢量，N； $\tau$ 为切应力张量。

c) 组分输运方程

$\frac{\partial \left(\rho c_f\right)}{\partial t}+\text{div}\left(\rho u{c}_f\right)=\text{div}\left[D_f\text{grad}\left(\rho c_f\right)\right]$ (3)

式中： $c_f$ 为组分 $f$ 的体积浓度； $D_f$ 为组分 $f$ 的扩散系数。

d) 标准 $\kappa -\epsilon$ 湍流模型。湍动能方程、湍动能耗散率方程、湍流粘性系数及速度梯度引起的湍动能方程分别为

$\frac{\partial }{\partial t}\left(\rho k\right)+\nabla \left(\rho uk\right)=\left(\mu +\frac{{\mu }_{\text{t}}}{{\sigma }_{\kappa }}\right){\nabla }^{2}k+G_k-\rho \epsilon +S_k$ (4)

$\frac{\partial }{\partial t}\left(\rho \epsilon \right)+\nabla \left(\rho u\epsilon \right)=\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right){\nabla }^2\epsilon +\epsilon \left(C_{1\epsilon }{G}_k-C_{2\epsilon }\rho \text{ε}\right)/\kappa +S_{\epsilon }$ (5)

${\mu }_t=C_{\mu }\rho k^2/\epsilon$ (6)

$G_k={\mu }_t\frac{\partial u_i}{\partial u_j}\left(\frac{\partial u_i}{\partial u_j}+\frac{\partial u_j}{\partial u_i}\right)$ (7)

式中： $u$ 是气体速度矢量； $\mu$ 表示流体的分子动力粘度； ${\mu }_t$ 表示湍流粘性系数； $k$ 是湍流动能； $\epsilon$ 表示湍流动能耗散率； $G_k$ 表示湍动能的生成项； $C_{\mu },C_{1\epsilon },C_{2\epsilon }$ 均为湍流常数 $C_{\mu }=0.09,C_{1\epsilon }=1.44,C_{2\epsilon }=1.92$ ； ${\sigma }_k,{\sigma }_{\epsilon }$ 均为湍流普朗特数， ${\sigma }_k=1.0,{\sigma }_{\epsilon }=1.3$ ； $S_k,S_{\epsilon }$ 均为用户定义源项； $u_i,u_j$ 表示 $i,j$ 方向的速度分量。

2.3. 边界条件设置

采用ANSYS Fluent软件对地下综合管廊燃气舱空间内的管道泄漏过程进行CFD仿真模拟。具体设置见表1与表2。

Table 1. Boundary conditions of numerical simulation

表1. 数值模拟边界条件

Table 2. Numerical simulation parameter settings

表2. 数值模拟参数设置

续表

3. 模型验证

为确保数值模型的可靠性，本文采用方自虎等人[12]的实验数据对模型结果进行对比验证，以评估所建模型的可靠性和准确性。验证实验采用尺寸为长10 m、宽0.15 m、高0.15 m的矩形舱体作为实验对象，实验中设定天然气释放流速为每分钟4 L，泄漏口尺寸为10 mm × 10 mm。通过实验结果与模拟结果的对比分析，确认模型的精度和有效性。

图2为No. 16号气体传感器采集的CO₂体积分数的模拟结果与试验数据的对比情况，从中可以观察到气体体积分数随扩散时间的变化趋势。总体来看，模拟结果与试验数据具有较高的一致性，气体体积分数在扩散初期快速增加，随后逐渐趋于稳定。在浓度值方面，仿真结果与试验数据的浓度值平均相对误差为19.26%。在303 s时，最小误差为6.55%。

综合上述分析，本研究提出的CFD仿真气体扩散模型满足在模拟综合管廊场景下气体泄漏与扩散过程中的适用性和精度要求。

Figure 2. No. 16 comparison of monitoring points results

图2. No. 16监测点结果对比图

4. 网格敏感性分析

在计算流体力学(CFD)模拟中，网格无关性分析是确保模拟结果可靠性和准确性的关键步骤。其主要目的是验证计算结果是否对网格划分的细化程度敏感，即通过逐步细化网格并比较计算结果的变化，判断是否已达到足够的网格分辨率，以捕捉流场中的主要物理特性。网格过于粗糙可能导致模拟结果偏离真实值，而过度细化则会显著增加计算成本但对精度提升有限。通过网格无关性分析，研究者能够在精度与计算效率之间实现平衡，从而确保所得结果具有科学性和工程应用价值。

4.1. 网格划分

网格划分方式直接影响计算效率和精度。因此，在泄漏位置为Z = 30 m，管道压力为1 MPa，入口风速为0.5 m/s，出口设置为表压排气，泄漏直径为100 mm的情况下，进行网格无关性研究。对网格数分别为406,972、431,615、468,668、518,180和587,419的五种结构化网格方案进行了评估。

4.2. 结果分析

图3为5种网格方案下，Z = 34 m、Z = 66 m、Z = 98 m监测点甲烷摩尔浓度变化图。

(a) z = 34 (b) z = 66

(c) z = 98

Figure 3. Table of change of methane molar concentration with different grid number

图3. 不同网格数甲烷摩尔浓度变化表

在图3(a)中，不同网格数量方案下，甲烷摩尔浓度的增长趋势在扩散初期(约0~15秒)存在一定的差异，特别是网格数量较少的方案(如318,180和387,419)与更高网格数量的方案相比，表现出略微不同的增长速率。然而，当扩散时间超过15秒后，各方案的浓度变化曲线逐渐趋于一致，且在最终阶段(40秒以后)，不同网格的浓度值基本重合。这表明网格数量对初期局部流场特性的捕捉存在一定影响，但对长期浓度场的描述效果相对一致。

在图3(b)中，不同网格方案在初期的曲线差异与图3(a)类似，但网格数量更高(如518,180和587,419)的方案在扩散早期(0~10秒)的结果更加接近实际的稳定趋势，说明较高网格数量能够更准确地捕捉到流体扩散的初始特性。随着扩散时间的增加，各曲线在15秒后逐步趋于一致，且最终结果的偏差可以忽略不计。这种趋势与图3(a)一致，进一步验证了网格划分的适用性。

图3(c)展示了离泄漏点更远的位置处的浓度变化。在初期扩散阶段(0~15秒)，与较低监测点相比，各网格方案之间的差异更小，这可能是由于高处气体扩散较均匀，流场特性对网格敏感性降低的原因。即使在网格数量较低的情况下(如406,972和431,615)，结果也能较早地趋于与高网格方案一致。这表明，在距离泄漏点较远的位置处，当前网格划分方案已能充分满足模拟的需求，进一步细化网格对结果影响较小。

综上，在网格数量达到468,668时，各监测点的甲烷浓度变化曲线在全时间范围内基本一致，表明模拟结果已经达到网格无关性。尽管在扩散初期，不同网格数量方案的模拟结果存在一定差异，但随着扩散时间增加，各方案的结果逐渐重合。此外，从空间位置的角度看，网格无关性的验证结果在不同位置处表现出一致性。因此，在后续的模拟中选择468,668个网格数。

4.3. 模拟结果

按照上述条件，通过ANSYS Fluent软件对管廊燃气舱内天然气管道的气体泄漏进行CFD仿真模拟，并对Fluent的仿真结果进行后处理，以结果云图形式呈现，更加直观地展示气体扩散的动态过程及其关键特征。具体仿真参数设置如下：进气入口风速为3 m/s，天然气管道压力为3 MPa，甲烷泄漏位置设定为Z = 90 m。

仿真结果通过云图展示了燃气舱内甲烷浓度的分布和扩散趋势(如图4所示)。从云图中可以清晰地观察到气体泄漏后的扩散特性，包括泄漏点附近的高浓度区域逐渐向四周扩散的过程，以及受气流驱动形成的浓度梯度分布。通过云图可直观反映气体泄漏扩散的关键特征及变化趋势，为后续分析燃气舱内气体泄漏风险和优化通风设计提供了科学依据。

Figure 4. Numerical simulation results

图4. 数值仿真结果

5. 结论

本文基于南宁市新邕路综合管廊燃气舱中的天然气泄漏工况，建立了三维CFD模型，对泄漏过程中甲烷的扩散行为进行了数值模拟，并通过实验数据验证了模型的可靠性。在此基础上，对网格无关性进行了深入研究，得出以下结论：

(1) 通过与文献中实验数据的对比，本文所建立的CFD模型在预测气体浓度随时间变化的趋势及量值方面均具有较高的精度。实验验证结果表明，仿真数据与实验数据的相对误差小于10%，可满足复杂场景中气体扩散模拟的准确性需求。

(2) 在网格数量达到468,668时，各监测点的甲烷浓度变化曲线在全时间范围内均趋于一致，表明模型已经达到网格无关性。进一步增加网格数量虽能略微提升结果精度，但对计算效率造成显著影响。因此，选择468,668个网格作为最优网格划分方案，在保证精度的同时有效降低计算成本。

(3) 研究表明，扩散初期(0~15秒)甲烷浓度的变化对网格数量较为敏感，尤其在靠近泄漏源的区域。然而，随着扩散时间的增加，不同网格方案的模拟结果逐渐趋于一致。高位监测点(如Z = 98 m)的结果表明，远离泄漏源的区域中，气体扩散更为均匀，对网格划分的敏感性较低，进一步验证了网格划分的适用性。

基金项目

国家自然科学基金资助项目(42077435, 42377171)；2022年南宁市科学研究与技术开发计划项目(ZC20223242)。

NOTES

^*第一作者。E-mail: 1067747181@qq.com

^#通讯作者。

参考文献