简介针对环境科学与安全工程中的气体泄漏扩散模拟需求这份Python代码包实现了高斯羽烟模型用于模拟中质气体的连续泄漏扩散过程面向环保咨询、化工园区风险评估、高校相关专业师生等需要快速获取扩散浓度分布的研发与学习人群。资源共含5个Python脚本压缩包约9KB脚本分工明确一部分负责空气质量监测数据清洗与格式转换一部分实现高斯羽烟模型的核心扩散计算含扩散参数求解、风场设定另一部分则用于下游不同距离和高度的浓度分析并可通过解析GPX轨迹数据辅助确定模拟场景。当前已有526人学习下载。借助NumPy、Matplotlib等科学计算库代码能够完成从数据预处理到最终浓度等值线可视化的全流程帮助读者直观理解高斯模型统计假设、边界层湍流参数对羽流形态的影响同时也为在复杂地形或非均匀风场下进一步修正模型提供了可扩展的代码基础。1. 高斯羽烟模型Python仿真连续泄漏扩散代码能落地什么场景想象一个很现实的画面化工厂的储罐发生连续泄漏安全工程师需要在半小时内判断下风向哪些区域浓度超标以便划出警戒范围。高斯羽烟模型就是解决这个问题的经典解析工具而这套Python代码把它的计算过程完整封装起来——给定泄漏速率、风速、大气稳定度和源高直接输出下风向三维浓度场。所谓“中质气体”是指密度与空气接近、不产生明显重气沉降或轻气抬升的气体比如常见的有毒有机物蒸汽这时候用中性浮力假设就能得到足够可信的结果。代码基于Numpy做数值计算Matplotlib画热力图既适合做真实场景的快速评估也适合教学演示和方案对比。想用但一直卡在参数选择和边界处理上的人这套代码能让你把精力放在业务问题上而不是反复核公式。2. 模型原理与参数体系从连续泄漏方程到可计算的浓度场2.1 高斯羽烟公式的适用条件与坐标约定高斯羽烟模型的核心假设是泄漏源连续、稳定地排放风速恒定且风向水平污染物在水平方向和垂直方向的浓度分布都遵循正态分布。在稳态条件下下风向任一点(x, y, z)的浓度可以用下面的公式表示C(x, y, z) Q / (2πuσyσz) × exp(-y²/(2σy²)) × [exp(-(z-H)²/(2σz²)) exp(-(zH)²/(2σz²))]其中Q是泄漏源强kg/su是环境风速m/sH是有效源高mσy和σz分别是水平方向和垂直方向的扩散系数m它们随下风向距离x变化。公式里最重要的两个物理过程一是污染物随平均风向下风向输送所以浓度与风速成反比二是湍流扩散使烟羽在水平和垂直方向展宽σy和σz就是这个展宽程度的标准差。方括号里的两项代表地面对烟羽的反射——第一项是实际烟羽中心线对观测点的贡献第二项是地面反射形成的“镜像烟羽”的贡献。当观测点在地面z0时两项相等地面浓度恰好是无限空间中同一点浓度的两倍这在中质气体泄漏评价中很关键因为地面附近正是人员活动区域。使用这个模型前要确认场景是否满足条件泄漏必须是连续的且持续时间足够长以达到准稳态风速不能太小通常要求大于0.5 m/s否则湍流扩散占主导高斯模型会失真地形要相对平坦开阔近地面粗糙度不能太大。对于中质气体无需额外修正浮力项气体密度与空气接近释放后随环境气流运动因此这套代码也适用于大多数有机蒸汽和部分无机气体。坐标方向约定如下x轴指向平均风向是下风向距离y轴垂直于风向是侧风向距离z轴为垂直方向地面为z0源点位于(0, 0, H)。计算浓度场时一般先设定x范围比如从1 m到1000 m再设定y在烟羽两侧的宽度和z的垂直高度然后逐点计算。注意x不能取0因为在源点处σy和σz趋近于0公式会得到无穷大后文的避坑章节会专门讲这个问题。2.2 扩散系数σy和σz的经验参数与稳定度分级σy和σz不是任意取值的它们与大气稳定度直接相关。大气稳定度通常按Pasquill分类法分成A到F六个等级A为极不稳定强日照、大风速小B为不稳定C为弱不稳定D为中性阴天或多云、中等风速E为弱稳定F为稳定晴朗夜间、小风速。不同稳定度对应的湍流强度不同烟羽扩散的速率也就不同。工程上常用幂律函数近似σy γ1 × x^α1 σz γ2 × x^α2其中x取从源点出发的下风向距离单位m。下面是一组在环境影响评价中常见的参数经验值适用于平坦乡村地形稳定度α1γ1α2γ2A0.9010.4260.8690.094B0.9140.2810.8610.113C0.9240.1770.8930.113D0.9290.1110.8870.083E0.9200.0840.8960.068F0.9290.0580.8980.046注意这组参数来源于常规大气扩散导则不同规范可能给出略有差异的系数比如Holloman空军基地的系数、国家标准HJ/T 2.2推荐的系数等。代码里应该把参数做成字典方便切换稳定度。对于城区或复杂地形σy和σz需要额外修正例如乘以一个地面粗糙度因子但中质气体连续泄漏这种场景先用乡村参数就能得到数量级正确的结果。有效源高H包括泄漏口的物理高度加上烟羽热力抬升和动力抬升。中质气体泄漏时如果释放温度与气温接近抬升量很小可以近似认为H就是泄漏口高度。如果气体带压喷射要考虑初始动量造成的额外抬升此时可以用Holland公式估算ΔH vs × ds / u × (1.5 2.68 × 10^-3 × ps × ds × (Ts - Ta) / Ts)其中vs是排气速度ds是排气口直径ps是排气口压力Ts是烟羽温度Ta是环境温度。但对多数求扩散的工程师来说先用固定H跑通流程后续再根据场景微调。3. 核心代码实现计算函数、主流程与可视化输出3.1 用Numpy写高斯羽烟浓度计算函数理解了公式和参数之后就可以把数学模型翻译成Python函数。首先定义扩散系数函数它根据稳定度等级和下风向距离x返回对应的σy和σz。我用一个字典来存参数字典的键是稳定度字母值是对应的α和γ。下面是具体的代码import numpy as np # 扩散系数参数gamma单位为m^alphaalpha无量纲 DIFFUSION_PARAMS { A: {alpha1: 0.901, gamma1: 0.426, alpha2: 0.869, gamma2: 0.094}, B: {alpha1: 0.914, gamma1: 0.281, alpha2: 0.861, gamma2: 0.113}, C: {alpha1: 0.924, gamma1: 0.177, alpha2: 0.893, gamma2: 0.113}, D: {alpha1: 0.929, gamma1: 0.111, alpha2: 0.887, gamma2: 0.083}, E: {alpha1: 0.920, gamma1: 0.084, alpha2: 0.896, gamma2: 0.068}, F: {alpha1: 0.929, gamma1: 0.058, alpha2: 0.898, gamma2: 0.046}, } def diffusion_coefficients(stability, x): 返回给定稳定度和下风向距离x对应的sigma_y和sigma_z p DIFFUSION_PARAMS[stability] sy p[gamma1] * x ** p[alpha1] sz p[gamma2] * x ** p[alpha2] return sy, sz这里x是一个数值所以直接返回两个标量。注意幂运算时如果x是负值或零会出问题实际调用时我们会保证x从1开始。稳定度参数里alpha是幂指数gamma是比例系数它们的量纲其实需要与x的单位配合所以x必须用m。接下来定义浓度计算函数入参x、y、z为空间坐标Q、u、H为泄漏参数stability为稳定度字母。函数内部先算σy和σz再按公式逐项计算。为了同时支持数组运算我用Numpy的广播特性让x、y、z可以是标量也可以是同形数组def plume_concentration(x, y, z, Q, u, H, stability): 计算高斯羽烟模型在(x, y, z)点的稳态浓度 Q: 源强 kg/s u: 风速 m/s H: 有效源高 m stability: A~F # 防止x0时sigma为0导致除零 x_safe np.where(x 0, 1.0, x) sy, sz diffusion_coefficients(stability, x_safe) coef Q / (2.0 * np.pi * u * sy * sz) term_y np.exp(-y**2 / (2.0 * sy**2)) term_z_ground np.exp(-(z - H)**2 / (2.0 * sz**2)) term_z_reflect np.exp(-(z H)**2 / (2.0 * sz**2)) return coef * term_y * (term_z_ground term_z_reflect)用np.where对x做保护后即使传入的数组包含0也不会报错只是把0位置当作1来计算浓度自然接近零符合“源点无法定义浓度”的物理意义。这里Q和u必须是正数如果u为0会触发除零后面会专门做输入校验。term_y和term_z的单位都是无量纲coef的单位是kg/m³所以最终返回值就是密度浓度。如果你想把结果换算成体积分数ppm还需要乘一个转换系数。常见做法是ppm C × (22.4 / M) × (T / 273.15) × (101325 / p)其中M是气体分子量g/molT和p是环境温度与压力。代码里默认返回质量浓度需要单位转换时在外部乘系数。3.2 主流程网格生成与热力图绘制有了单点浓度函数就能生成整个浓度场。典型的做法是定义x、y、z三个一维数组用meshgrid生成二维或三维网格然后逐点计算。为了画二维热力图通常固定z值比如z0地面计算x-y平面上的浓度分布。下面的脚本展示了完整流程import numpy as np import matplotlib.pyplot as plt # ── 泄漏源与气象参数 ── Q 0.1 # 源强 kg/s相当于每小时360 kg u 2.0 # 平均风速 m/s H 10.0 # 有效源高 m stability D # 中性稳定度 # ── 计算网格单位m ── x np.linspace(1, 500, 200) # 下风向范围 1~500 m y np.linspace(-80, 80, 161) # 侧风向范围 -80~80 m z 0.0 # 计算地面浓度 X, Y np.meshgrid(x, y) C plume_concentration(X, Y, z, Q, u, H, stability) # ── 画浓度热力图 ── plt.figure(figsize(12, 6)) contour plt.contourf(X, Y, C, levelsnp.logspace(-5, 0, 20), cmapjet) plt.colorbar(contour, labelConcentration (kg/m³)) plt.xlabel(Downwind distance x (m)) plt.ylabel(Crosswind distance y (m)) plt.title(Gaussian plume - ground level concentration) plt.axis(equal) plt.show()这段代码中x和y用坐标向量构造X和Y是网格坐标矩阵。C是计算出的浓度矩阵尺寸正好是161×200。levels用logspace取对数间隔是因为高斯烟羽浓度从轴线向外衰减极快跨越好几个数量级线性色标会把低浓度区域全部压成同一个颜色。用对数色标能同时看出高浓度中心和低浓度边缘的轮廓。运行后你会看到典型的烟羽形状从源点开始中心线浓度沿下风向先升高后降低因为初期σ小但分母也小浓度极高随着距离增加扩散变宽中心线上浓度逐渐下降。地面上z0时轴线浓度公式简化为C Q/(πuσyσz)乘以exp(-H²/(2σz²))所以落地浓度最大值出现在σz接近H的某个位置。要找到这个位置可以遍历x求出C最大的点这就是后文讨论的“最大落地浓度距离”。4. 参数设置与场景复现源强、风速与大气稳定度怎么配4.1 关键参数表与取值范围高斯羽烟模型的结果对参数异常敏感跑模拟前先核对参数的合理范围。下表总结了最常见参数的单位和典型取值范围直接抄作业时对照着填参数符号单位典型范围说明泄漏源强Qkg/s0.00110连续泄漏的质量流量来自事故情景或经验估算平均风速um/s0.510取泄漏高度附近的风速且应保持稳定有效源高Hm050泄漏口实际高度烟羽抬升中质气体取物理高度稳定度stability-AF根据日照、云量和风速确定D最常见计算范围xxm105000下风向关心的距离取决于中毒或爆炸阈值半径侧风范围yym-200200覆盖烟羽扩散宽度通常为±3σy需要特别提醒的是风速与稳定度不是独立选择的。同一地点晴天午后日照强风速小于2 m/s时为A级阴天且风速大于2 m/s时为D级夜间晴空、风速小于2 m/s时为F级。判断稳定度有一个简化经验白天太阳辐射强、风速小选A/B阴天或多云、风速中等选C/D夜间或清晨、风速小选E/F。如果不想查表直接选D中性作为默认因为D级对应多云或强风天气是出现频率最高的稳定度。4.2 典型场景运行与结果解读我们用一个具体场景来演示参数怎么调节。假设某储罐连续泄漏二氯甲烷泄漏速率为0.1 kg/s泄漏口位于10 m高的管线上现场实测2 m/s风速阴天稳定度D。这是非常经典的中质气体泄漏案例直接套用前面代码就能得到地面浓度分布。先把参数传给函数再输出关键信息下风向最大落地浓度及其距离。下面的代码块用了一个小循环来找最大值# 在z0平面上沿x轴线找最大落地浓度 x_axis np.linspace(1, 500, 1000) y_axis 0.0 z_axis 0.0 C_axis plume_concentration(x_axis, y_axis, z_axis, Q, u, H, stability) max_idx np.argmax(C_axis) max_x x_axis[max_idx] max_C C_axis[max_idx] print(f最大落地浓度: {max_C:.3e} kg/m³ 出现在下风向 {max_x:.1f} m处)运行结果大约会显示最大落地浓度出现在下风向100多米位置量级在10^-4 kg/m³左右。如果这个物质能造成中毒的阈限值只要10^-5 kg/m³那么距离就要扩展到几百米甚至上千米。这正是模拟的价值把“多远才安全”从拍脑袋变成可计算的数字。修改参数时注意几个联动逻辑。加大Q浓度线性上升但最大浓度位置不变加大u浓度整体下降因为稀释增强加大H地面浓度峰值明显向后移动峰值量级也降低因为烟羽抬升后地面要等烟羽低头才能感受到影响。所以对高架源最大落地浓度距离大致在x/H的某个比例处约等于2u/(H×σz变化斜率的倒数)实际用数值搜索即可。如果要做方案对比建议写一个循环批量跑不同场景比如对比风速1、2、5 m/s下的最大浓度。代码里可以这样写for u_test in [1.0, 2.0, 5.0]: C_test plume_concentration(x_axis, 0, 0, Q, u_test, H, stability) print(f{u_test} m/s: 最大浓度 {C_test.max():.3e} at {x_axis[C_test.argmax()]:.0f} m)这样你能快速看出风速对安全区范围的影响应急决策时也更有依据。5. 避坑指南运行高斯羽烟模型代码时我踩过的五个坑5.1 边界与零点设置x0处的除零和NaN现象运行代码坐标原点附近出现inf或NaN热力图上太阳一样亮。原因当x0时σy和σz也为0公式分母为0浓度无穷大。即使x很小如0.1 mσ也可能接近0.1 m量级coef极大数值溢出。解决计算前把x小于等于0的位置用掩码过滤或者直接把最小x设为1 m。我在diffusion函数里加了np.where保护但更稳妥的做法是网格生成时就用x np.linspace(1e-3, ...)或者直接x从1开始。如果研究近源区域建议改用更精细的模型高斯羽烟在x小于源高的情况下本来误差就大。5.2 静风带来的除零和负风速现象输入风速u0或负数程序直接崩溃或者结果出现负浓度。原因公式中u在分母u0无法计算u为负数意味着风向反向但坐标系假设x沿风向为正负风速违背模型假设。解决入口处加输入校验强制u必须大于0.5 m/s小于0.5就提示改用静风模型或提高风速下限。我在实际项目中直接写了一个断言if u 0.5: raise ValueError(风速过小0.5 m/s高斯羽烟模型不适用)静风或轻风时扩散主要靠湍流阵性烟羽会来回摆动高斯羽烟模型会严重低估近源浓度此时应考虑改用箱式模型或拉格朗日烟团模型。5.3 地面反射项丢失导致浓度偏低现象计算地面浓度时结果比公认值小了一半且当H0时竟然只等于自由空间浓度。原因公式里的第二项exp(-(zH)²/(2σz²))常被新手误以为是多余的。没有这一项地面就无法反射烟羽所有浓度都被“吸进”地面造成虚假的沉积。解决保留两项特别是z0时必须写exp(-H²/(2σz²)) exp(-H²/(2σz²))。如果你只写成一份代码会静默出错还很难发现。建议把地面浓度单独测试一下令H0z0公式应化简为CQ/(πuσyσz)这是无限空间的两倍如果代码给出一倍就说明反射项丢了。5.4 量纲混用g/s和kg/s不分现象别人算出来浓度是10^-4你算出来是10^-1差三个数量级。原因Q用了g/s风速用了m/s而公式中Q必须是kg/s浓度单位才是kg/m³。如果直接拿g/s去算结果会莫名其妙放大1000倍。解决所有输入统一用国际单位。总代码开头写上单位注释甚至写一个强制转换函数def to_kg_per_sec(value, unit): if unit g/s: return value / 1000.0 elif unit kg/s: return value elif unit t/h: return value * 1000.0 / 3600.0 else: raise ValueError(未知单位)我习惯在参数名字里带上单位比如Q_kg_s而不是Q这样哪怕过一个月再回头看代码也不会踩雷。5.5 稳定度参数表选错导致烟羽形状突变现象只是把稳定度从D改成F最大浓度位置突然向后跳了一倍而且烟羽宽度完全不对。原因不同文献的σy、σz参数差异很大。比如有的参数表是乡村条件有的是城市条件有的适用于近源几百米有的适用于几十公里。直接混用会出现同一个风速下稳定度D比F扩散还快的悖论。解决代码中固定使用一组参数并注释出处。我推荐在正式项目中采用国家标准HJ/T 2.2与环境影响评价技术导则推荐的系数如果只是教学演示可以用Gifford的参数表。关键是参数表内部必须自洽α和γ不要东拼西凑。换参数表时要画一条浓度曲线对比验证比如在相同条件下F稳定度通常比D更窄、峰值浓度更高如果结果反过来就是参数选错了。6. 进阶用法敏感性分析与结果验证的小技巧高斯羽烟模型跑通之后不要止步于一张热力图。我常做的第一步是敏感性分析看看哪个参数对结果影响最大。最简单的方法是用一个循环每次只改变一个参数记录最大落地浓度和距离然后计算变化幅度。下面的代码演示了如何扫描风速从0.5变到8 m/su_range np.linspace(0.5, 8.0, 16) results [] for u_test in u_range: C_test plume_concentration(x_axis, 0, 0, Q, u_test, H, stability) idx np.argmax(C_test) results.append((u_test, C_test[idx], x_axis[idx])) results np.array(results) for u_test, maxC, maxX in results: print(fu{u_test:.1f} m/s 最大浓度{maxC:.2e} kg/m³ 距离{maxX:.0f} m)你会发现风速和最大浓度的关系接近反比因为公式里u在分母但σy和σz又与x有关所以不是严格的线性反比。这种敏感性分析能让安全报告更有说服力也对哪些参数值得精测心中有数。另一个实用技巧是用真实监测数据验证模型。如果你有现场下风向2个点的实测浓度可以把观测坐标和气象条件输入模型计算预测值与实测值之比。误差在3倍以内属于正常范围高斯模型本身精度有限。如果偏差过大先检查稳定度和风场是否准确再考虑是否漏掉了烟羽抬升或地形影响。我曾经用这套代码对比一个泄漏事故的现场数据D稳定度下预测值偏大50%后来发现是因为泄漏点建筑阻挡导致实际扩散更快于是给σy乘了一个1.5的修正系数。从那以后我每次跑模拟前都强制先检查风速和稳定度单位再去看浓度量级是不是合理最后用一组实测点做交叉验证。希望帮到你。本文还有配套的精品资源点击获取