简介气体泄漏扩散模拟是环境风险评估与应急决策的核心技术其本质是通过数学模型描述污染物在风场中的输运与稀释规律。高斯烟羽模型因其计算高效、参数物理意义明确成为工业事故后果模拟和数学建模竞赛中的主流方法。该模型基于湍流扩散理论结合大气稳定度、源强、风速等参数可快速生成下风向浓度场分布为疏散范围划定和应急预案提供定量依据。在核电站泄漏场景中借助Python实现浓度场计算与等值线可视化能够有效衔接理论公式与工程实践支撑从泄漏源项分析到应急分区决策的完整链条。本文围绕此类问题的建模流程详解模型选型、参数处理与结果验证帮助读者掌握气体扩散模拟的落地路径。1. 核电站泄漏后的放射性气体扩散从应急浓度场到数学建模竞赛的完整解法核电站一旦发生泄漏第一时间要回答的问题不是“泄漏了多少”而是“下风向哪个位置、在什么高度、什么时间会出现多高的放射性浓度”。这个问题的答案直接决定疏散范围、隔离区划和应急决策也是各类数学建模竞赛里最经典的“源项—气象—扩散—浓度场”四段式题目。本文围绕“43367核电站泄漏后放射性气体浓度分布规律和气体扩散模型”这个赛题方向完整拆解从机理模型选型、参数表设计、Python 计算、可视化到结果验证的落地路径。这套方法不仅适用于华为杯、全国大学生数学建模竞赛这类赛事也适用于辐射环境影响评价、化工厂有毒气体泄漏后果模拟等工业场景。我会把“高斯烟羽模型怎么建模”“扩散参数查哪张表”“等值线图怎么画才不画错”“为什么你的浓度场总是被评委质疑”这几个关键问题一次讲透。新手照步骤能跑通老手能直接拿去改参数换场景。2. 气体扩散模型的数学框架高斯烟羽与烟团模型怎么选2.1 高斯模型的适用前提为什么核电站泄漏首选它而不是 CFD放射性气体泄漏的扩散模拟工业界和竞赛里最常用的是高斯型扩散模型而不是 CFD计算流体力学。原因很简单CFD 要解 N-S 方程、要画网格、要算湍流闭合一个泄漏事故场景跑一次要数小时而高斯模型在分钟级就能给出满足应急决策精度的浓度场。高斯烟羽模型的核心假设是连续泄漏的污染气体在平稳均匀的风场中浓度在横向y和垂直向z都服从高斯分布向下风向x方向按风速输送并逐渐稀释。浓度场的解析表达式为C(x, y, z) Q / (2π · u · σ_y · σ_z) · exp(-y² / (2σ_y²)) · [exp(-(z-H)² / (2σ_z²)) exp(-(zH)² / (2σ_z²))]其中 Q 是源强Bq/su 是平均风速m/sH 是有效排放高度mσ_y 和 σ_z 是水平和垂直扩散参数m随下风向距离 x 和大气稳定度变化。这个公式里的“地面反射项”exp(-(zH)²/(2σ_z²)) 模拟地面全反射——气体分子撞到地面不掉进土里而是弹回大气。那什么时候用烟羽模型、什么时候用烟团模型判断标准是泄漏方式和气象条件。持续泄漏连续几小时甚至几天用烟羽模型短时泄漏、爆炸式释放、或者风速风向变化剧烈要用烟团模型——把泄漏源在时间轴上切成一系列瞬时释放的“烟团”每个烟团单独扩散再在空间点上叠加。华为杯、国赛题里给“事故发生在某日某时持续泄漏 2 小时”这种条件就属于典型的可分段处理的准连续释放。2.2 扩散参数 σ_y 和 σ_z这是整个模型里最容易被扣分的环节σ_y、σ_z 不是常数它们随下风向距离 x 和大气稳定度等级A~F 六级变化。常见的参数化方案有 Pasquill-Gifford 曲线、Briggs 公式和国家标准推荐的幂函数形式。Briggs 模式因为形式简洁、竞赛中可以直接查表带公式用得最多。Briggs 扩散参数开阔乡村地表稳定度等级σ_y (m)σ_z (m)A0.22x(10.0001x)^(-1/2)0.20xB0.16x(10.0001x)^(-1/2)0.12xC0.11x(10.0001x)^(-1/2)0.08x(10.0002x)^(-1/2)D0.08x(10.0001x)^(-1/2)0.06x(10.0015x)^(-1/2)E0.06x(10.0001x)^(-1/2)0.03x(10.0003x)^(-1)F0.04x(10.0001x)^(-1/2)0.016x(10.0003x)^(-1)注意这张表的适用距离范围是 0.1 到 10 公里。超过 10 公里σ_z 的增长会变缓甚至饱和因为大气混合层高度限制了垂直方向的进一步发展。很多参赛队伍把公式硬套到 50 公里外画出来的浓度等值线在 20 公里处出现明显的“喇叭口收缩”——评委一眼就觉得你的模型没有考虑混合层顶盖效应。稳定度等级怎么定国标 GB/T 3840-91 里有一套太阳高度角、云量、风速查表法。但竞赛题通常不会给这么细的气象数据常见的做法是直接按题给条件判断晴天白天且风速小于 2 m/s 定为 B 级阴天或多云天定为 D 级晴天夜间且风速小定为 F 级。如果题里给的是“夜间、风速 1.5 m/s”那就是 F 级对应的 σ_z 公式分母里带 (10.0003x)下风向 5 公里处垂直扩散只有几十米浓度衰减极慢——这种场景下放射性烟流会贴着地面传输很远应急疏散范围要比白天大得多。3. 从泄漏源到浓度场Python 实现、参数表和边界条件3.1 最小可运行代码单点连续泄漏浓度场计算拿到赛题第一步是把“物理描述”翻成“可计算条文”。下面这份代码是我最常用的一套最小骨架输入源强、风速、稳定度、泄漏高度输出一个二维浓度场z 取地面高度。import numpy as np import pandas as pd import matplotlib.pyplot as plt def briggs_sigma(x, stability): Briggs 扩散参数x: 下风向距离数组(m) if stability A: sy 0.22 * x / np.sqrt(1 0.0001 * x) sz 0.20 * x elif stability B: sy 0.16 * x / np.sqrt(1 0.0001 * x) sz 0.12 * x elif stability C: sy 0.11 * x / np.sqrt(1 0.0001 * x) sz 0.08 * x / np.sqrt(1 0.0002 * x) elif stability D: sy 0.08 * x / np.sqrt(1 0.0001 * x) sz 0.06 * x / np.sqrt(1 0.0015 * x) elif stability E: sy 0.06 * x / np.sqrt(1 0.0001 * x) sz 0.03 * x / (1 0.0003 * x) elif stability F: sy 0.04 * x / np.sqrt(1 0.0001 * x) sz 0.016 * x / (1 0.0003 * x) return sy, sz def gaussian_plume(x, y, z, Q, u, H, stability): 高斯烟羽模型返回浓度场 C (Bq/m^3)x沿风向 X, Y np.meshgrid(x, y) sy, sz briggs_sigma(X, stability) # 地面反射项: 源高 H 与镜像源 -H 叠加 term_z np.exp(-(z - H)**2 / (2 * sz**2)) np.exp(-(z H)**2 / (2 * sz**2)) C Q / (2 * np.pi * u * sy * sz) * np.exp(-Y**2 / (2 * sy**2)) * term_z return C # 情景参数源强 3.7e16 Bq/s约 100 万居里/小时量级 # 风速 4 m/s有效释放高度 60 mD 级中性稳定度 Q 3.7e16 u 4.0 H 60.0 stability D x np.linspace(100, 30000, 300) # 下风向 0.1~30 km y np.linspace(-5000, 5000, 200) # 横向 ±5 km z 0.0 # 地面浓度 C gaussian_plume(x, y, z, Q, u, H, stability) # 浓度超过 1e3 Bq/m^3 的区域面积km^2 threshold 1e3 area np.sum(C threshold) * (x[1]-x[0]) * (y[1]-y[0]) / 1e6 print(f浓度 {threshold} Bq/m^3 的覆盖面积: {area:.2f} km^2)这份代码的核心逻辑有三个第一briggs_sigma函数把稳定度等级翻译成随距离变化的扩散参数第二浓度公式里的term_z同时算了实际源和地面镜像源的贡献保留这一项才能算出地面浓度第三网格只覆盖下风向 30 公里、横向 ±5 公里的矩形区域够辐射应急前期判断用。3.2 参数表设计把题目给的数据映射到模型输入赛题不会直接给你Q3.7e16这么整的参数给你的通常是“堆芯放射性活度总量 1.2e18 Bq泄漏率按总量的 0.1% 每小时估算”“事故持续释放 4 小时后切断”。这种情况下源强的换算要分两步先算每小时泄漏量再均匀分摊到每秒钟作为 QBq/s。下面是核电站泄漏场景最常用的参数映射表输入项赛题常见给法模型取值方法源强 Q“堆芯总活度 1.2e18 Bq泄漏率 0.1%/h”Q 1.2e18 × 0.001 / 3600 ≈ 3.33e11 Bq/s有效释放高度 H“安全壳破裂高度 40 m热释放抬升”抬升高度用 ΔH 1.6·u_f^(-2/3)·(Q_h)^(1/3)或题给直接值风速 u“10 m 高度风速 4 m/s”换算到释放高度u u_10 × ln(H/0.1) / ln(10/0.1)稳定度“夜间少云风速 2 m/s”查稳定度等级表得 F少数取 E混合层高度不直接给默认 1000 mD 级夜间可压到 300~500 m有效释放高度是另一个常见的翻车点。有热源时烟气会因热浮力和出口动量向上抬升不能直接拿安全壳破口高度当 H。最简单的抬升公式是 Holland 公式ΔH (1.5·w·d/u) × (1 ΔT/T_s)其中 w 是出口流速d 是破口直径ΔT 是烟气与环境温差。实际竞赛里如果题目没给破口流速和温度就直接用 H 破口标高并在灵敏度分析里说明 H 对浓度结果的影响幅度——这个处理比硬套 Holland 公式更稳妥。3.3 坐标旋转把计算坐标系对齐真实风向模型公式里的 x 永远是沿风向的轴。但赛题给的地图是正北朝向泄漏源在 (x0, y0)监测点的经纬度或平面坐标是 (xi, yi)。直接拿地图坐标代入公式就错了——必须先把所有监测点坐标旋转到沿风向下风向一维、侧风向一维的坐标系。def rotate_coords(xs, ys, wind_angle_deg): 将地图坐标(正北为y轴正东为x轴)旋转到沿风向坐标系。 wind_angle_deg: 风向方位角0北风90东风即风从哪个方向来。 返回 (x_along_x, y_cross)x_along 沿风向y_cross 垂直风向。 # 风向方位角转弧度风向北偏东 wind_angle_deg 度 theta np.radians(90 - wind_angle_deg) x_along (xs - xs[0]) * np.cos(theta) (ys - ys[0]) * np.sin(theta) y_cross -(xs - xs[0]) * np.sin(theta) (ys - ys[0]) * np.cos(theta) # 只保留下风向点 mask x_along 0 return x_along[mask], y_cross[mask], mask # 示例泄漏源在(5000, 5000)风向为正西风270°从西往东吹 source np.array([5000.0, 5000.0]) monitor_x np.array([5500.0, 9000.0, 3000.0]) monitor_y np.array([5200.0, 4800.0, 5000.0]) x_along, y_cross, mask rotate_coords( monitor_x - source[0], monitor_y - source[1], wind_angle_deg270 )这个代码里最容易出错的是“风向方位角”的定义气象上“北风”指风从北面吹来对应风向方位角 0°风往南吹而坐标系旋转需要的是“风吹去的方向”。所以代码里用90 - wind_angle_deg做了换算。如果你的监测点在源的上风向x_along 0模型算出的浓度是 0——浓度不为 0 但极其小这个在应急判断里可以直接忽略但在论文里最好写一句“上风向浓度可忽略”体现你处理过这个边界。4. 扩散规律可视化与定量解读等值线图怎么画才对4.1 烟羽轴向浓度衰减先画一维曲线验证量级直接扑向二维等值线图之前建议先画一维轴向浓度衰减曲线。这个步骤能暴露 90% 的参数错误——比如 σ_z 公式输入了负的 x、Q 量级少了 3 个零、稳定度等级输反了。# 沿烟羽中心线y0, z0的浓度 x_axis np.linspace(100, 30000, 500) sy_a, sz_a briggs_sigma(x_axis, stability) C_axis Q / (np.pi * u * sy_a * sz_a) * np.exp(-(0 - H)**2 / (2 * sz_a**2)) # 基础量级检查近源 1 km 处浓度应在 1e5 ~ 1e7 Bq/m^3 范围 print(fx1km C{C_axis[np.argmin(abs(x_axis-1000))]:.2e} Bq/m^3) print(fx10km C{C_axis[np.argmin(abs(x_axis-10000))]:.2e} Bq/m^3) plt.figure(figsize(10, 5)) plt.semilogy(x_axis/1000, C_axis, b-, labelfstability{stability}) plt.axhline(1e3, colorr, ls--, lw1, label1e3 Bq/m^3) plt.xlabel(下风向距离 (km)) plt.ylabel(地面浓度 (Bq/m^3)) plt.legend() plt.grid(True, alpha0.3) plt.show()这段检查代码的价值在于D 级稳定度、4 m/s 风速下1 公里中心线浓度如果在 1e6 量级附近说明参数正常如果出现 1e20 或者 1e-8基本就是单位换算错了Bq/s 和 Bq/h 没转干净或者 Q 用错了单位。对数坐标下浓度曲线应该是一条整体平顺下降的线任何“鼓起”“凹陷”都说明 σ_y、σ_z 在某个距离段出现了明显不连续——多数情况是np.sqrt(1 0.0001 * x)里的 x 用了公里单位但公式要求米。4.2 地面浓度等值线风向坐标系画完再转回地图坐标系等值线图是竞赛论文里最直观的成果图也是最容易被挑刺的地方。画图流程分三步先在高斯坐标系里画出以源为起点、沿风向拉长的浓度场再把坐标旋回正北地图坐标系最后叠加地形、河流、居民区底图。# 在沿风向坐标系中计算浓度场 x_mesh np.linspace(100, 15000, 200) y_mesh np.linspace(-3000, 3000, 200) C_mesh gaussian_plume(x_mesh, y_mesh, 0, Q, u, H, stability) plt.figure(figsize(12, 4)) levels [1e1, 1e2, 1e3, 1e4, 1e5, 1e6] cs plt.contourf(x_mesh/1000, y_mesh/1000, C_mesh, levelslevels, cmaphot, normLogNorm()) plt.colorbar(cs, label浓度 (Bq/m^3)) plt.xlabel(下风向距离 (km)) plt.ylabel(侧风向距离 (km)) plt.title(地面浓度等值线风向坐标系) plt.show()画完等值线要把坐标轴刻度从米换算成公里等值线层级用levels[1e1, 1e2, ...]这种十进制间隔方便竞审判读时估算半致死区通常取 1e5 Bq/m^3 对应的范围、限入区1e3 Bq/m^3和警戒区1e1 Bq/m^3。千万别把contour和contourf混用后不设置色标范围——默认的 10 级彩色会把 1e1 和 1e5 压在同一色阶里图面看起来整片都是红色毫无区分度。用LogNorm()做颜色归一化是关键一步。5. 核电站泄漏建模的避坑清单稳定度、源项和坐标系的 5 个血泪教训5.1 稳定度等级选错导致疏散范围差 5 倍现象同一场泄漏白天和夜间两种方案的结果差了一整个数量级下风向 10 公里的地面浓度白天算出来是 10 Bq/m^3夜间算出来是 800 Bq/m^3。原因白天太阳辐射强大气不稳定垂直混合快污染物能快速扩散到高空地面浓度低夜间地表辐射冷却形成逆温层垂直扩散被压制放射性气体被困在近地面层地面浓度极高。很多队伍忽略了“夜间”这个题眼统一按 D 级中性条件算。解决严格按题给时间、云量和风速查稳定度表。夜间少云风速小于 2 m/s 直接定为 F 级如果题目给的是阴雨天才可以用 D 级。竞赛答案里疏散半径的计算结果会因为一个稳定度等级的差别变动 3~5 倍这是评委非常看重的区分点。5.2 高斯分布侧风向浓度出现负值坐标现象代码输出的浓度场在源的上风向出现非零浓度或者侧风向超过 10 公里的区域仍有大片高浓度区。原因rotate_coords里坐标偏移量计算少了源坐标分量或者风向角换算时三角函数符号搞反。更隐蔽的是 y 的网格范围取的太窄±2 km而 F 级稳定度下侧风向扩散极慢2 km 外浓度仍然显著截断效应让等值线在边界处密集堆积。解决先用正北风向wind_angle0做一次自检确认浓度场中心线沿 y 轴正向分布再用正东风向wind_angle90验证是否沿 x 轴正向分布。侧风向网格至少取到 σ_y 的 6 倍——F 级在 30 公里处 σ_y 约 1.2 公里侧风网格至少要 ±8 公里。5.3 源强 Q 的量纲混乱Bq/s 与 Bq/h 的灾难现象算出的 1 公里地面浓度高达 1e15 Bq/m^3明显违反物理直觉。原因赛题给的泄漏率如果是“每小时 3.7e16 Bq”你直接喂给模型当 Q相当于把每秒泄漏量放大了 3600 倍。这是新手最常见的黑匣子问题——公式里每个符号都有明确单位但代入时没人核查。解决写代码时在参数表上方留一行单位注释# Q: Bq/s不是 Bq/h计算中途打印一次 Q 的值和量级最终结果核验时用“总释放活度守恒”检查对全空间积分浓度 × 风速 × 面积应该等于 Q。5.4 远距离扩散的混合层顶盖效应未处理现象30 公里外的地面浓度比实测或预期高 2~3 倍且等值线越拉越长像一根细尾巴。原因高斯烟羽模型的 σ_z 随距离无限增长但实际大气有混合层顶盖通常在 300~1500 m污染物碰到顶盖后无法继续向上扩散只能水平摊开。近源 5 公里内顶盖效应不明显超过 10 公里后误差显著。解决远距离模拟用“多次反射镜像法”替代单个镜像项——在 σ_z 超过混合层高度 H_mix 之前把 z 方向高斯项改写成无限镜像源叠加或者直接用分段公式σ_z ≤ H_mix/1.5 时用单镜像超过之后改用均匀混合假设。竞赛中如果不做这么细的修正至少要在论文里写一句“本研究忽略混合层顶盖效应结果对 10 km 以远浓度存在高估约 30%~50%”——坦诚限界比硬拗更显专业。5.5 有效释放高度 H 被误用为地理高度现象把安全壳高度 40 米直接当 H 代入公式结果地面浓度比用 H60 米的方案高了近一个量级。原因泄漏气体带热热浮力使烟羽中心线上抬。40 米处的破口因为加热抬升实际“有效高度”可达 80~120 米。H 越高地面浓度越低但影响范围更远。解决如果题目给了热释放速率 Q_hMW用 Briggs 抬升公式估算 ΔH没给就从保守角度取小值并在敏感性分析中把 H 从 0.5 倍到 2 倍区间扫一遍画出浓度最值和面积随 H 的变化曲线。这个图放在论文里几乎必加分。6. 浓度场结果的敏感性分析与竞赛论文的验证闭环模型跑通只是开始要说服评委“你的模型可靠、结果可信”必须做敏感性分析和验证闭环。我自己的习惯是每套浓度场都必须过四道检查源强守恒、监测点倒算、稳定度敏感性、疏散范围对比。第一道守恒检查最简单对任意垂直剖面做通量积分——以源为圆心、半径 R 的圆柱面上总通量应等于 Q。代码里用C * u * dA沿下风向 10 公里垂直截面求和和 Q 的偏差超过 15% 说明网格分辨率不够或边界截断太狠。第二道是监测点倒算。赛题通常会给出若干个环境监测点在事故后某时刻的读数把监测点坐标代入模型算出预测浓度与实测值对比计算偏差因子模拟值/实测值。因子在 0.3~3 之间说明模型可用偏差超过 10 倍优先怀疑源强估算而不是扩散参数——源强反演本来就是数学建模竞赛里一个高频考点用正演模型配最小二乘反解 Q。第三道敏感性分析参数矩阵至少包含风速0.5、1、2、4、8 m/s和稳定度D、F两个维度。从我的实证经验看风速对地面浓度的影响近似反比风速翻倍浓度减半而稳定度从 D 到 F 会让 10 公里中心线浓度提升 3.5~6 倍。这个量级认知能帮你在赛题要求“给出不同气象情景下的应急响应分级”时快速答出有依据的范围。最后一关也是我踩过的那道坎把浓度阈值转成地图上的实际疏散区域。输出不能停留在“浓度等值线图”要落到“安全区、控制区、限制区”的边界坐标用表格列出每个分区的面积和典型居民点是否包含在内。国赛和华为杯这类赛题在阅卷时很看重“从模型结果到决策建议”这条链的闭环只给图不给决策建议等于把后半题丢了。建议的进阶扩展方向是加一个简单的烟团示踪动画——把泄漏过程切成 5 分钟一个烟团每个烟团按当前时刻的风速和稳定度独立扩散叠加显示地面浓度随时间演化。这个思路能把静态等值线变成动态响应过程也符合题目里“泄漏发生后 12 小时内浓度分布规律”这类时间演化要求。希望这套从模型选型、参数表设计、Python 实现到验证闭环的路径能帮你把核电站泄漏扩散这个赛题做成竞赛里的加分项。做这类题最舒服的状态不是把所有方法都堆上去而是把高斯模型吃透、参数说明白、边界交代清楚最后用一张等值线图和一张疏散分区表把故事讲完整。希望帮到你。本文还有配套的精品资源点击获取