简介本资源是一套基于Python实现的高斯羽烟模型气体扩散模拟代码面向环境科学、安全工程及应急响应领域的初学者与实践者用于定量模拟中质气体在大气中的连续泄漏扩散过程。压缩包共5个文件全部为Python脚本.py总大小仅9KB轻量易用其中gpm_1.py与gpm_2.py为核心计算模块实现高斯羽烟浓度场求解downstream_look.py支持下游关键位置浓度追踪分析convert-aqms.py提供空气质量监测数据格式转换能力gpx-parser.py可解析GPS轨迹以辅助风场参数设定。已有525人学习下载代码结构清晰、注释充分配套参数文件如123Y-2、G2已内嵌典型工况设置开箱即可运行并快速验证泄漏影响范围是开展风险评估、预案推演与教学演示的实用工具集。1. 项目概述从“烟”到“羽”的工业安全模拟在化工、能源、消防应急乃至城市公共安全领域一个核心且紧迫的问题是当有毒或易燃气体发生泄漏时它究竟会扩散到哪里浓度有多高多久会影响到周边的居民区或敏感设施十年前回答这些问题可能依赖昂贵的计算流体力学CFD软件和大型服务器集群门槛极高。但现在任何一个掌握基础Python编程的安全工程师、环境评估师甚至学生都能借助经典的高斯羽烟模型快速构建一个可靠的气体扩散模拟工具。这正是我们这个项目的核心价值将复杂的工业安全与风险评估转化为一段清晰、可复现的Python代码。高斯羽烟模型常与高斯烟团模型一同被提及是大气扩散模拟的基石。简单理解你可以把泄漏源想象成一个持续冒烟的烟囱羽烟或瞬间爆炸释放的烟团烟团。模型通过一系列基于大量实测数据归纳出的数学公式描述这些“烟”在下风向的浓度分布。它虽然做了许多理想化假设如均匀风速、稳定气象条件但其计算速度快、参数物理意义明确、结果易于解读在事故快速预警、安全距离划定、环境影响评价等场景中有着不可替代的实用价值。本次我们聚焦于“连续泄漏”也就是模拟像管道小孔持续逸出气体这样的情景这正是高斯羽烟模型的典型应用。通过这篇分享你将获得一个完整的、从中质气体如天然气、氨气、氯气等密度接近空气扩散行为相对典型泄漏模拟到可视化分析的全流程代码实现。我会带你一步步拆解公式、编写代码、处理边界条件并分享我在实际风险评估项目中积累的调试技巧和避坑指南。无论你是想完成课程设计、开展科研预研还是为实际工作构建一个快速分析工具这篇文章都能提供可直接“抄作业”的解决方案。2. 模型核心原理与关键假设拆解在动手写代码之前我们必须吃透模型背后的物理图像和它的“能力边界”。高斯模型之所以高效正源于其一系列精心设计的假设理解这些假设你才能知道何时该用它何时结果可能不靠谱。2.1 高斯羽烟模型的基本方程高斯羽烟模型的核心公式描述的是在下风向任意一点(x, y, z)处的气体浓度C。对于地面点源连续泄漏其经典形式如下C(x, y, z) (Q / (2π * u * σy * σz)) * exp(-y² / (2 * σy²)) * [exp(-(z - H)² / (2 * σz²)) exp(-(z H)² / (2 * σz²))]这个公式看起来复杂但我们可以像拆解机器一样把它分解开Q 源强。这是模拟的“起点”单位通常是kg/s或g/s。它表示泄漏源每秒释放多少质量的物质。准确获取或估算Q是整个模拟可信度的关键通常基于泄漏孔径、压力、流体性质等计算得出。u 风速。模型假设风速是恒定且均匀的方向沿x轴正方向。这是模型最大的理想化之一实际中风是湍流、变化的。σy和σz 水平和垂直方向的扩散参数。这是高斯模型的“灵魂”。它们不是常数而是随着下风向距离x的增加而增大的函数即σy(x)和σz(x)。它们代表了烟羽在传播过程中因为大气湍流而不断向四周散开的程度。σ值越大烟羽越“胖”浓度分布越分散。H 有效源高。对于地面泄漏H通常为0。如果泄漏点有抬升如烟囱H就是物理高度加上烟羽的抬升高度。指数项exp(...) 这部分描述了浓度在y横向和z垂直方向上的高斯正态分布形态。中心轴线上浓度最高向两边呈钟形曲线衰减。公式中最后两个exp项相加是考虑了地面的反射作用相当于在真实源下方-H处放置了一个镜像源使得地面附近浓度更高。注意 公式中[exp(-(z - H)²/(2σz²)) exp(-(z H)²/(2σz²))]这一项正是处理地面反射的标准方法。它保证了在地面z0处垂直方向的浓度梯度符合物理规律。如果你的模拟不关心垂直分布只关心地面浓度可以令z0公式会简化。2.2 扩散参数σy, σz的选取模型的“经验库”扩散参数σy和σz如何确定它们无法从第一性原理推导而是基于海量的野外示踪实验数据拟合出来的经验公式。最常用的体系是帕斯奎尔-吉福德Pasquill-Gifford, P-G曲线及其对应的数学拟合公式。P-G体系首先根据风速、太阳辐射和云量等气象观测将大气稳定度划分为A到F六个等级A极不稳定 晴朗夏日午后湍流强烈扩散极快。B不稳定、C弱不稳定 常见的白天条件。D中性 阴天或大风夜晚最常用也最“平均”的假设。E弱稳定、F稳定 晴朗少风的夜晚大气层结稳定扩散能力弱容易导致近地面高浓度积聚是最危险的气象条件之一。每个稳定度等级都对应一组σy(x)和σz(x)的曲线或公式。在代码中我们通常采用 Briggs 等人提出的城市/乡村条件下的分段幂律公式来近似这些曲线例如对于乡村条件的D类中性稳定度σy 0.08 * x / sqrt(1 0.0001 * x)σz 0.06 * x / sqrt(1 0.0015 * x)这里的x是下风向距离。选择与你的模拟场景城市/乡村和气象条件最匹配的稳定度等级是影响结果准确性的最重要决策之一。2.3 模型的适用性与局限性高斯模型是一个“宝藏工具”但绝非“万能钥匙”。它的核心假设决定了其应用边界稳态条件 风速、风向、泄漏源强在模拟期间恒定不变。这显然与瞬息万变的现实有差距因此它更适用于短时、定常的泄漏情景分析。均匀流场 忽略地形、建筑物引起的复杂风场变化。在平坦开阔地带模拟效果较好在山谷或建筑群中误差会增大。无沉降、无化学反应 模型假设泄漏物质是惰性的在扩散过程中不与空气发生化学反应也不因重力沉降或雨水冲刷而损失。这对于许多中质气体是合理的但对于重气如液化石油气初始阶段或反应性气体如氯气遇水则需要更复杂的修正模型。适用于一定范围 通常适用于下风向几公里到几十公里的范围。太近如百米内的复杂近场扩散或太远百公里以上的输送模型误差较大。实操心得 在实际风险评估中我通常用高斯模型做“第一轮筛查”或“保守估计”。例如在未知具体气象条件时我会选择最不利的F类稳定度扩散最慢和较低的风速进行计算这样得到的地面最大浓度和影响范围是偏保守的用于划定安全红线非常有效。如果需要更精确的厂区内部扩散CFD模拟才是更好的选择但耗时和成本也呈指数级增长。3. Python代码实现从公式到可视化理解了原理我们就可以开始搭建代码骨架了。我们将采用模块化设计使代码清晰、易复用、易扩展。3.1 环境准备与依赖库我们只需要基础的数值计算和绘图库。建议使用 Anaconda 创建一个新的虚拟环境。# 创建并激活环境可选但推荐 conda create -n gas-diffusion python3.9 conda activate gas-diffusion # 安装核心库 pip install numpy matplotlibNumPy 处理数组运算的核心所有网格计算和公式实现都依赖它。Matplotlib 用于绘制浓度等高线图、三维曲面图等直观展示扩散结果。3.2 核心函数编写高斯浓度计算我们将核心计算封装成一个函数输入位置坐标和气象参数输出浓度值。import numpy as np def gaussian_plume_conc(x, y, z, Q, u, H, stability_classD, ruralTrue): 计算给定点的高斯羽烟模型浓度。 参数: x, y, z : float or array 接收点坐标 (m)。x为下风向距离y为横风向距离z为高度。 Q : float 源强 (kg/s)。 u : float 风速 (m/s)。 H : float 有效源高 (m)。地面泄漏设为0。 stability_class : str 大气稳定度等级A到F。 rural : bool True为乡村条件False为城市条件。影响扩散参数公式。 返回: C : float or array 浓度 (kg/m^3)。通常转换为 mg/m^3 或 ppm 以便于分析。 # 1. 定义扩散参数σy, σz与下风向距离x的函数关系以Briggs公式为例 # 这里仅以乡村D类为例实际需要实现一个完整的查找或计算函数 def get_sigma(x, stability, rural): # 这是一个简化示例你需要根据P-G表或Briggs公式扩展此函数 # 为所有稳定度等级和城乡条件提供系数a, b, c, d # 公式形式通常为σ a * x / (1 b * x)**c 或 σ d * x**e if stability D and rural: sigma_y 0.08 * x / np.sqrt(1 0.0001 * x) sigma_z 0.06 * x / np.sqrt(1 0.0015 * x) elif stability F and rural: # 稳定条件扩散更慢 sigma_y 0.04 * x / np.sqrt(1 0.0001 * x) sigma_z 0.016 * x / np.sqrt(1 0.0003 * x) else: # 默认回退到中性条件 sigma_y 0.08 * x / np.sqrt(1 0.0001 * x) sigma_z 0.06 * x / np.sqrt(1 0.0015 * x) return sigma_y, sigma_z # 确保x不为零或负值对于源点附近需特殊处理 x np.maximum(x, 1.0) # 避免除以零 # 2. 获取当前x处的扩散参数 sigma_y, sigma_z get_sigma(x, stability_class, rural) # 3. 应用高斯公式 # 第一部分系数项 coefficient Q / (2 * np.pi * u * sigma_y * sigma_z) # 第二部分横向y方向高斯分布 exp_y np.exp(-0.5 * (y / sigma_y) ** 2) # 第三部分垂直z方向高斯分布包含地面反射 exp_z1 np.exp(-0.5 * ((z - H) / sigma_z) ** 2) exp_z2 np.exp(-0.5 * ((z H) / sigma_z) ** 2) # 4. 计算浓度 C coefficient * exp_y * (exp_z1 exp_z2) return C # 单位: kg/m^3关键点解析get_sigma函数是代码的“数据心脏”。上面的示例极度简化。在实际应用中你需要根据选定的 Briggs 系数表或其它权威来源实现一个完整的、支持所有 P-G 稳定度等级的查找函数。这通常是一个包含大量if-elif或字典查找的代码块。公式中的np.maximum(x, 1.0)是一个小技巧。因为扩散参数公式在x0附近可能不适用或导致计算溢出将其限制为一个小的正值如1米可以保证计算的稳定性且对稍远距离的结果影响微乎其微。浓度单位是kg/m^3。对于有毒气体我们更常用mg/m^3或ppm体积分数。需要进行单位转换C_mg_per_m3 C * 1e6从 kg/m^3 到 mg/m^3转到ppm需要知道气体的摩尔质量M(g/mol) 和标准状态下的摩尔体积约22.4 L/mol公式为C_ppm (C_mg_per_m3 * 22.4) / M。这个转换在后续结果分析中非常重要。3.3 构建模拟区域与网格计算为了可视化我们需要在整个关心的地理区域网格上计算浓度。def create_simulation_grid(x_range, y_range, z0, dx10, dy10): 创建用于模拟的二维平面网格。 参数: x_range : tuple 下风向距离范围例如 (0, 2000) 单位米。 y_range : tuple 横风向距离范围例如 (-500, 500)。 z : float 计算平面高度通常为0地面浓度。 dx, dy : float 网格分辨率 (m)。越小越精细计算量越大。 返回: X, Y : 2D arrays 网格点的坐标矩阵。 x np.arange(x_range[0], x_range[1] dx, dx) y np.arange(y_range[1], y_range[0] - dy, -dy) # 注意顺序保证绘图方向正确 X, Y np.meshgrid(x, y) return X, Y def calculate_concentration_grid(X, Y, Z, Q, u, H, stabilityD, ruralTrue): 在整个网格上计算浓度。 参数: X, Y, Z : 2D arrays 网格坐标。 ... 其他参数同 gaussian_plume_conc ... 返回: C_grid : 2D array 每个网格点上的浓度 (kg/m^3)。 # 初始化一个与网格同形状的零数组 C_grid np.zeros_like(X) # 获取网格形状 rows, cols X.shape # 遍历每个网格点向量化操作避免低效循环 # 这里我们利用NumPy的广播和数组运算直接对整个数组进行计算 # 注意我们的 gaussian_plume_conc 需要能接受数组输入 # 修改原函数使其内部运算支持数组上面的示例已基本支持 C_grid gaussian_plume_conc(X, Y, Z, Q, u, H, stability, rural) return C_grid实操心得 网格分辨率dx, dy的选择需要权衡。对于初步筛查50米甚至100米的分辨率可能就够了。但如果要精确找到最大浓度点或绘制光滑的等高线可能需要10米或更高的分辨率。一个技巧是可以先粗算定位高浓度区域再在该区域进行局部网格加密计算以节省总计算时间。3.4 结果可视化让数据“说话”可视化是分析模拟结果的关键。我们将创建地面浓度等高线图和沿中心轴线的浓度剖面图。import matplotlib.pyplot as plt def plot_ground_level_contour(X, Y, C_grid, Q, thresholdNone): 绘制地面浓度等高线图。 参数: X, Y : 2D arrays 网格坐标。 C_grid : 2D array 浓度网格数据 (kg/m^3)。 Q : float 源强用于标题。 threshold : float, optional 关注的最低浓度阈值会绘制一条加粗的等值线。 C_mg_per_m3 C_grid * 1e6 # 转换为 mg/m^3 plt.figure(figsize(10, 8)) # 绘制等高线 contour_levels np.logspace(np.log10(C_mg_per_m3.min()1e-6), np.log10(C_mg_per_m3.max()), 15) cp plt.contourf(X, Y, C_mg_per_m3, levelscontour_levels, cmapviridis, alpha0.75) plt.colorbar(cp, label浓度 (mg/m³)) # 绘制特定阈值等值线如IDLH浓度 if threshold is not None: threshold_mg threshold * 1e6 # 假设输入的threshold单位是 kg/m^3 cs plt.contour(X, Y, C_mg_per_m3, levels[threshold_mg], colorsred, linewidths3, linestyles--) plt.clabel(cs, fmtf{threshold_mg:.1f} mg/m³, inlineTrue, fontsize10) # 标记泄漏源位置 plt.scatter(0, 0, colorblack, s100, marker*, label泄漏源, zorder5) plt.axhline(y0, colorgray, linestyle-, linewidth0.5) # 中心线 plt.xlabel(下风向距离 (m)) plt.ylabel(横风向距离 (m)) plt.title(f高斯羽烟模型模拟 - 地面浓度分布 (Q{Q} kg/s)) plt.legend() plt.grid(True, alpha0.3) plt.axis(equal) # 保证x和y轴比例相同图形不变形 plt.tight_layout() plt.show() def plot_centerline_profile(X, C_grid): 绘制下风向中心轴线y0的浓度变化曲线。 参数: X : 2D array 网格x坐标。 C_grid : 2D array 浓度网格数据。 # 找到y0对应的行索引假设网格对称中心行 center_row_idx C_grid.shape[0] // 2 C_centerline C_grid[center_row_idx, :] * 1e6 # 转换为 mg/m^3 x_centerline X[center_row_idx, :] plt.figure(figsize(10, 5)) plt.plot(x_centerline, C_centerline, b-, linewidth2, label中心轴线浓度) plt.xlabel(下风向距离 (m)) plt.ylabel(浓度 (mg/m³)) plt.title(下风向中心轴线浓度剖面) plt.grid(True, alpha0.3) plt.legend() # 标记最大浓度点 max_conc_idx np.argmax(C_centerline) max_conc C_centerline[max_conc_idx] max_dist x_centerline[max_conc_idx] plt.scatter(max_dist, max_conc, colorred, s80, zorder5) plt.annotate(f最大值: {max_conc:.2f} mg/m³\n距离: {max_dist:.0f} m, xy(max_dist, max_conc), xytext(max_dist*1.1, max_conc*0.8), arrowpropsdict(arrowstyle-, colorred)) plt.tight_layout() plt.show()4. 完整模拟流程与参数化分析现在我们将所有模块组合起来运行一个完整的模拟案例并进行参数敏感性分析。4.1 一个完整的模拟案例假设我们模拟一个化工厂的氨气NH₃储罐发生小孔连续泄漏。氨气分子量约为17 g/mol属于中质气体。# 模拟参数设置 Q 0.1 # 源强0.1 kg/s一个中等规模的泄漏 u 2.0 # 风速2 m/s微风条件 H 0 # 地面泄漏 stability D # 中性稳定度 rural True # 乡村环境 # 创建模拟网格下风向2公里横风向左右各500米 x_range (0, 2000) y_range (-500, 500) X, Y create_simulation_grid(x_range, y_range, z0, dx20, dy20) # 计算浓度场 print(正在计算浓度场...) C_grid calculate_concentration_grid(X, Y, 0, Q, u, H, stability, rural) print(计算完成。) # 可视化 # 设置一个关注阈值例如氨气的立即危害生命健康浓度(IDLH)约为 270 ppm。 # 需要转换为质量浓度。粗略估算270 ppm NH₃ ≈ 270 * (17/24.45) ≈ 188 mg/m³ (24.45为25°C时摩尔体积L/mol) # 对应 kg/m^3 为 1.88e-4 threshold_conc 1.88e-4 # kg/m^3 plot_ground_level_contour(X, Y, C_grid, Q, thresholdthreshold_conc) plot_centerline_profile(X, C_grid) # 输出关键信息 C_max C_grid.max() C_max_mg C_max * 1e6 print(f最大地面浓度: {C_max_mg:.2f} mg/m³) # 找到最大浓度点的位置 max_idx np.unravel_index(np.argmax(C_grid), C_grid.shape) x_max X[max_idx] y_max Y[max_idx] print(f最大浓度点位置: (x{x_max:.0f} m, y{y_max:.0f} m))运行这段代码你将得到两张图一张是彩色的地面浓度分布等高线图红色虚线标出了IDLH浓度阈值范围另一张是沿下风向中心线的浓度变化曲线清晰地展示了浓度如何先上升后下降并标出了峰值位置和数值。4.2 参数敏感性分析风速与稳定度的影响高斯模型的结果对输入参数非常敏感。我们通过简单的循环来量化这种影响。# 分析不同风速和稳定度的影响 wind_speeds [1.0, 2.0, 5.0] # m/s stability_classes [B, D, F] # 不稳定中性稳定 # 固定其他参数 Q_fixed 0.1 H_fixed 0 x_eval 1000 # 评估下风向1000米处的浓度 y_eval 0 z_eval 0 results [] for u in wind_speeds: for stab in stability_classes: C gaussian_plume_conc(x_eval, y_eval, z_eval, Q_fixed, u, H_fixed, stab, ruralTrue) results.append({ 风速 (m/s): u, 稳定度: stab, 浓度 (mg/m³): C * 1e6 }) # 用表格展示结果 import pandas as pd df_results pd.DataFrame(results) print(\n不同条件下下风向1000米处中心轴线浓度) print(df_results.pivot(index稳定度, columns风速 (m/s), values浓度 (mg/m³)).round(2))你会得到一个类似下面的表格稳定度1.0 m/s2.0 m/s5.0 m/sB (不稳定)较低值较低值最低值D (中性)中等值我们案例的值低值F (稳定)最高值高值中等值分析结论一目了然风速的影响风速越大稀释作用越强同一点位浓度越低。风速减半浓度可能近似翻倍。稳定度的影响在静稳天气F类大气垂直混合能力弱污染物不易向上扩散导致近地面浓度极高是最危险的气象条件。而不稳定条件B类有利于垂直扩散地面浓度较低。这个简单的分析突显了在风险评估中考虑“最坏情况”低风速、稳定层结的重要性。5. 常见问题、调试技巧与模型扩展在实际编码和应用中你肯定会遇到各种问题。这里分享一些我踩过的坑和解决思路。5.1 数值计算问题与稳定性处理在源点附近出现“尖峰”或异常高值原因 当x趋近于0时扩散参数σy和σz也趋近于0导致高斯公式中的分母趋于零计算溢出或产生不现实的极高浓度。解决 如前所述对x设置一个最小值如1.0米。物理上模型本身在近源处通常100倍泄漏口径就不适用。因此模拟结果在源点附近是无效的应忽略。浓度等高线图出现“锯齿”或不光滑原因 网格分辨率dx, dy太粗或者扩散参数公式在某个距离处有拐点。解决 提高网格分辨率。对于等高线绘制可以尝试在plt.contourf中使用更精细的levels或者先对C_grid进行轻微的平滑处理如使用scipy.ndimage.gaussian_filter但需谨慎避免过度失真。单位混淆导致结果差几个数量级典型错误 源强Q用了g/s但以为结果是kg/m^3或者风速u用了km/h未转换成m/s。检查清单 始终使用国际单位制SI。Q: kg/s,u: m/s,x,y,z: m输出C: kg/m^3。在最终呈现时再转换为mg/m^3或ppm。5.2 模型局限性及实用修正高斯模型是理想化的针对一些特定情况有经验性的修正方法有限混合层高度 在实际大气中逆温层会像一个“盖子”限制污染物的垂直扩散。当烟羽高度H 2.15σz接近混合层高度L时污染物会被限制在混合层内反复反射导致地面浓度增加。此时高斯公式需要引入一个包含无穷级数的修正项来模拟多次反射。对于距离较远的计算这是一个重要的修正。非地面源高架源 我们的代码中H参数就是用来处理这个的。对于烟囱H是物理高度加上烟羽的抬升高度。烟羽抬升计算本身就是一个复杂的子课题有霍兰德Holland等多种经验公式。重气扩散的初始阶段 像液化石油气LPG这种冷重气体泄漏初期会像“瀑布”一样下沉并形成重气云团其扩散行为完全不符合高斯模型。此时需要使用专门的重气模型如SLAB、DEGADIS等。一个粗略的修正方法是在重气主导阶段使用一个“虚拟源”概念即认为经过一段距离的初始塌陷和空气卷吸后重气云团转变成了中性气体再从那个虚拟源位置开始使用高斯模型。实操心得 对于应急响应速度往往比绝对精度更重要。一个在1分钟内给出保守估计结果的高斯模型比一个需要1小时设置和运行的精细CFD模型更有价值。因此理解模型的保守性偏向例如在稳定气象条件下预测的浓度更高并善用它是将其应用于实际安全决策的关键。5.3 代码优化与扩展方向向量化与性能 我们的代码已经利用了NumPy的广播机制性能对于中等网格如200x200足够了。如果网格非常大可以考虑使用numba的jit装饰器对核心计算函数进行即时编译可提升数十倍速度。集成更多功能多源叠加 如果有多个泄漏源总浓度场就是各个源单独产生的浓度场的线性叠加。只需循环计算每个源的贡献并累加即可。动态气象输入 可以编写一个包装函数读入随时间变化的风速、风向序列分段进行稳态模拟然后将结果在时间维度上拼接近似模拟非稳态过程。风险评估输出 将浓度场与毒性浓度阈值如IDLH、ERPG等比较自动计算并绘制出“影响区域”并计算该区域内的人口或资产暴露情况。构建图形用户界面GUI 使用PyQt或Streamlit可以快速构建一个参数输入面板和结果展示仪表盘让不熟悉代码的同事或客户也能使用这个工具。Streamlit尤其适合快速部署为Web应用。这个基于Python的高斯羽烟模型模拟项目就像为你打开了一扇通往大气扩散模拟和工业风险评估的大门。它轻量、透明、可定制让你能够深入理解每一个参数背后的物理意义。从这段代码出发你可以根据具体的项目需求不断打磨扩散参数表、添加修正算法、优化可视化最终将其打造成一个得心应手的专业分析工具。记住所有模型的输出都只是对现实的近似而工程师的价值就在于深刻理解这种近似的边界并在此基础上做出审慎可靠的判断。本文还有配套的精品资源点击获取