简介一份以核电站泄漏为背景的数学建模竞赛完整参赛文档面向需要研究气体扩散模型与应急影响评估的建模学习者。内容系统性覆盖连续源高斯烟羽模型、瞬时源一维与三维抛物型扩散模型、有限时间泄漏叠加模型以及有风条件下风向确定时的浓度预测方法结合福岛核泄漏案例展示空气扩散、食品与工业产品传播等实际影响评估思路并附有陕西师范大学2011年模拟赛题全文、模型假设、符号说明与MATLAB求解过程。包体为单份doc文档大小仅1.18MB便于完整阅读、对照修改和复用。该资源已有204人学习浏览包含完整的问题重述、模型建立与求解、模型改进等竞赛结构便于模仿论文组织与复用公式推导适合作为数学建模专题训练、核应急扩散模拟、环境风险评估或相关课程设计的参考资料。1. 一个核电站泄漏建模题为什么大家都在用高斯烟羽把核电站泄漏后放射性气体浓度分布规律和气体扩散模型研究这道数学建模竞赛题打开第一反应通常是补核物理知识但真正动手后会发现决定分数高低的不是辐射剂量学而是你对大气扩散模型的理解深度和使用精度。这类题目本质上要求你回答三件事——泄漏源有多强、气体往下风向怎么扩散、地面上的人会吸入多少。绝大多数队伍都会选择高斯烟羽模型作为基线它写得出解析解、参数含义清晰、画得出漂亮的等浓度线而且评委最容易验证你的中间步骤。这篇文章就按我给参赛队辅导时的完整思路来拆模型怎么选、公式怎么改、参数怎么定、哪些地方最容易翻车。2. 扩散模型选型烟羽、烟团与湍流模拟的适用边界2.1 三类模型的数学结构与成本对比处理大气污染物扩散常见做法有三条路高斯烟羽模型、高斯烟团模型以及基于拉格朗日粒子或计算流体力学的数值模拟。如果只做一次竞赛题高斯类是性价比最高的因为它把湍流扩散简化为浓度在横风向和垂直方向都服从高斯分布最后得到的是闭式解几十行代码就能跑出全场浓度。高斯烟羽模型适合描述连续泄漏。它的数学基础是稳态对流扩散方程假设源强稳定、风速恒定、湍流场均匀那么下风向任意一点的浓度只与该点坐标、源强、风速和扩散参数有关。烟团模型则是把一团放射性物质视为一个在风中漂移且不断膨胀的高斯分布云团适合瞬时泄漏。CFD 类模型能处理复杂地形和建筑绕流但网格划分、湍流闭合方案、边界条件每一项都能耗掉你大半比赛时间而且评委很难快速复现你的结果。所以我的建议非常明确竞赛论文里以高斯烟羽为主模型用烟团做瞬时泄漏对比把 CFD 留到改进方向里提两句不要真去跑一个非稳态湍流模拟。表格里三类模型的核心差异是这样。模型类型数学形式适用场景计算成本竞赛中的定位高斯烟羽解析解连续泄漏、稳态气象极低主模型高斯烟团解析解瞬时泄漏、泄漏后时间演化极低对比模型拉格朗日粒子随机游走积分复杂流场、非均匀地形中高改进方向CFD数值求解N-S方程建筑绕流、复杂地形高一般不碰2.2 麻烦的大气稳定度Pasquill 分级与扩散参数的取法高斯烟羽公式里最关键的输入是水平扩散参数 σy 和垂直扩散参数 σz它们不是常数而是下风向距离和处理大气稳定度的函数。大气稳定度决定了湍流混合的强弱不稳定大气里热力湍流旺盛烟羽很快被摊开稳定大气里湍流受到抑制烟羽又窄又贴地地面浓度反而更高。这是扩散建模里第一个反直觉结论。Pasquill 稳定度分级把大气分成 A 到 F 六类A 为极不稳定D 为中性F 为极稳定。判断依据是风速、太阳辐射强度和云量。竞赛题通常会直接给你「晴天、风速 3 m/s、中午」这类气象描述你要把它翻译成稳定度等级。白天按太阳辐射强弱分夜间按云量多少分这个查表过程我放在后文实操部分先记住一点同一道题A 类和 F 类算出的最大落地浓度能差一个量级选错等级等于白做。拿到稳定度等级后用 Briggs 公式计算扩散参数。Briggs 给出的是 σy、σz 随距离 x单位 km变化的经验式参数表里 A、B、C、D、E、F 各有一套系数。很多队伍直接照抄系数就完事却没有注意 x 的量纲导致浓度整体偏移。我习惯在代码里先把 x 统一换算成 km 再代入和 Briggs 公式的适用范围保持一致。2.3 连续泄漏还是瞬时泄漏先判断场景再选模型题目描述里「持续释放」「冷却剂持续流失」「释放持续 120 小时」这类措辞指向连续泄漏用烟羽模型如果是「安全壳压力骤升后瞬间破裂」「爆炸性抛射」这类措辞指向瞬时泄漏用烟团模型。还有一种折中情况释放持续了若干小时但你看的是泄漏后很短时间内几百米范围内的浓度这时烟团和烟羽差异明显。判断方法其实很朴素比较释放时长和烟气从源到关心点所需的时间。设释放时间为 T烟气到达距离 x 需要的时间是 x/u。如果 T 远大于 x/u可以近似为连续泄漏如果 T 与 x/u 同阶甚至更短必须按瞬时泄漏处理。竞赛里很多题目会把放射性气体释放描述成「事故后第 4 小时开始以稳定速率持续排放」这种情况直接当成连续源问题会被大大简化。还有一条经验当题目给出了泄漏持续时间、但不同核素释放速率可能变化时可以把时间轴切成几段每段当作一个稳态烟羽最后叠加浓度。这样做虽然有点粗糙但比用一个恒定源强全程计算要诚实得多评委也更容易认可你在工程近似上的判断。3. 把浓度分布写进程序高斯烟羽公式的完整落地与修正3.1 高斯烟羽公式拆解每个符号都对应一个物理过程连续点源高斯烟羽模型的浓度公式是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 分别是水平向和垂直向扩散参数x 是下风向距离y 是横风向距离z 是离地高度。我得提醒一点公式里的平方项决定了浓度分布的形状y 和 z 都以 m 为单位但 σy、σz 是用距离 x 算出来的单位也是 m。真实场景中如果出现 σz 为负或为零的输入程序会直接算出 NaN这种细节在数据预处理阶段就要拦下来。公式中的 [exp(-(z-H)²/2σz²) exp(-(zH)²/2σz²)] 包含了实源和地面反射镜像源两项。物理含义是烟羽向下扩散到地面时地面不会吸收放射性气体而是把它反射回大气数学上用镜像源叠加来模拟这个反射过程。缺失这一项靠近地面的网格浓度会明显偏低或偏高取决于你怎么处理边界条件。3.2 核素特有的两个修正衰变项和沉积扣除核电站泄漏的放射性气体和普通化工厂泄漏的毒气有个本质区别核素会衰变。设某核素的衰变常数为 λ单位 1/s气体从源传播到下风向距离 x 需要时间 t x/u沿程浓度应乘衰减因子 exp(-λx/u)。λ 用半衰期 T1/2 计算λ ln2/T1/2。以中等半衰期的放射性核素为例如果题目给出某核素半衰期为 8 天而关心的范围只有几十公里风速 3 m/s 下烟气到达最远点的时间不过几个小时衰变修正几乎可以忽略。但半衰期短的核素比如某些惰性气体同位素的半衰期只有几分钟这个因子能让几十公里外的浓度掉一到两个量级不能省。干沉积和湿沉积是第二类修正。干沉积指气溶胶粒子在地面、植被表面的沉降损失湿沉积指降雨清洗。竞赛题如果不给降雨信息一般不做湿沉积但干沉积可以用沉积速度 vd 近似处理把垂直项里的源项乘一个损耗因子或者在地面边界条件上引入沉积通量。多数获奖论文的做法是在模型改进部分讨论沉积影响主模型只在衰变项上做文章。这样既体现了物理完整性又不会被沉积参数的表征问题拖住。3.3 地面与混合层反射镜像源处理让浓度不虚高大气边界层上方存在一个混合层顶相当于第二块反射面。烟羽在混合层内上下反射多次数学上就是一组无限镜像源求和。实际计算中反射到三四阶已经足够收敛再多就是数值自嗨。具体做法是改写垂直项把 (z-H)、zH、z-H2zi、zH2zi 这些高度差全部加入高斯项zi 是混合层高度。如果简单地忽略混合层反射近地面的垂直扩散会被低估浓度计算结果在远距离处偏高等浓度线画出来会像一根被压扁的香肠形状明显失真。代码实现时可以单独写一个函数计算垂直项指数之和循环反射阶数从 0 到 3 或 4并设置一个阈值当新增项的贡献小于主项百分之一时停止。这种做法既保留了物理意义又把计算量控制在线性级别网格再密也跑得动。3.4 可运行的Python实现扫网格算浓度并找出最大落点下面是我在备赛时常给队伍用的最小实现。它接受源参数和气象参数返回三维浓度并自动找出最大地面浓度位置。import numpy as np def calc_sigma(x_km, pg_class): # x_km: 下风向距离, 单位 km # pg_class: A~F, Pasquill 稳定度等级 # 采用 Briggs 开放地形参数, 输入 x 单位需为 km x max(x_km, 0.01) # 防止 x0 时除零 # 参数来自 Briggs 经验公式, 使用时注意 x 单位为 km, sigma 单位为 m if pg_class A: sy 0.22 * x / np.sqrt(1 0.0001 * x) sz 0.20 * x elif pg_class B: sy 0.16 * x / np.sqrt(1 0.0001 * x) sz 0.12 * x elif pg_class C: sy 0.11 * x / np.sqrt(1 0.0001 * x) sz 0.08 * x / np.sqrt(1 0.0002 * x) elif pg_class D: sy 0.08 * x / np.sqrt(1 0.0001 * x) sz 0.06 * x / np.sqrt(1 0.0015 * x) elif pg_class E: sy 0.06 * x / np.sqrt(1 0.0001 * x) sz 0.03 * x / (1 0.0003 * x) else: # 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, pg_class, decay_lambda0.0, h_mix1000.0): # x: 下风向距离(m), y: 横风向距离(m), z: 高度(m) # Q: 源强(Bq/s), u: 风速(m/s), H: 有效源高(m) # decay_lambda: 衰变常数(1/s), h_mix: 混合层高度(m) if u 0 or H 0: raise ValueError(风速必须为正数, 源高不能为负) sy, sz calc_sigma(x / 1000.0, pg_class) if sy 0 or sz 0: raise ValueError(扩散参数必须为正数, 请检查稳定度等级与距离范围) # 衰变修正: 烟气从源到该点经历时间 t x / u decay_term 1.0 if decay_lambda 0 else np.exp(-decay_lambda * x / u) # 横风向高斯分布 lateral np.exp(-y**2 / (2.0 * sy**2)) # 垂直项: 实源 地面镜像 混合层顶反射镜像 vertical 0.0 # n 0 表示实源与地面镜像, n /-1, /-2 表示混合层顶反射 for n in range(-3, 4): z1 z - (H - 2 * n * h_mix) z2 z (H - 2 * n * h_mix) vertical np.exp(-z1**2 / (2.0 * sz**2)) vertical np.exp(-z2**2 / (2.0 * sz**2)) C (Q / (2.0 * np.pi * u * sy * sz)) * lateral * vertical * decay_term return C # 场景: 某核电站烟囱高度 120 m, 源强 1e9 Bq/s, 风速 4 m/s, 稳定度 D 类 params {Q: 1e9, u: 4.0, H: 120.0, pg: D, decay_lambda: 1e-6} # 计算地面浓度并找最大值 x_grid np.arange(100, 20000, 100) # 下风向 100m ~ 20km y_grid np.arange(-3000, 3001, 100) # 横风向 -3km ~ 3km max_c 0.0 max_pos (0, 0) for x in x_grid: for y in y_grid: c gaussian_plume(x, y, 0, params[Q], params[u], params[H], params[pg], params[decay_lambda]) if c max_c: max_c c max_pos (x, y) print(最大地面浓度: {:.4e} Bq/m^3.format(max_c)) print(出现位置: 下风向 {} m, 横风向 {} m.format(max_pos[0], max_pos[1]))这段代码的逻辑分三层。calc_sigma 函数把稳定度等级映射成 Briggs 扩散参数它是公式里唯一与气象有关的变量决定烟羽的胖瘦gaussian_plume 函数实现浓度计算内部包含了横向高斯分布、垂直镜像源求和、衰变修正最下面的双层循环扫描地面网格找到浓度最大值及其位置。垂直镜像循环从 -3 到 3实际就是在混合层顶上下各做了三次反射叠加程序里我在注释里说明了这一项的作用。参数上需要注意三处。第一calc_sigma 的输入 x 单位是 kmgaussian_plume 里传入的是 x/1000.0这样 Briggs 公式系数才不会出错第二decay_lambda 如果设为 0衰变项直接跳过一次 exp 计算避免不必要的数值开销第三双层循环在网格 200×61 时约运行 1.2 万次纯 Python 大约几秒如果题目要求反演或蒙特卡洛建议对循环做向量化用 numpy 的广播一次性算整个网格。4. 竞赛实操源项估算、气象参数与模型检验4.1 源项数据怎么拿从题目信息估算释放率竞赛题目一般不会直接给「源强1e9 Bq/s」这么舒服的数值常见的是给你堆芯总活度、安全壳泄漏率、释放窗口时间。这时候要自己搭估算链源强 Q 堆芯某核素总活度 × 释放份额 / 释放持续时间。举个例子题目假设某核素堆芯总活度为 1e17 Bq事故后安全壳完整性部分丧失泄漏份额约 0.01释放持续 24 h则 Q 1e17 × 0.01 / (24×3600) ≈ 1.16e10 Bq/s。这里面最容易错的点是把时间窗口和泄漏份额搞混。题目里的「持续泄漏」和「一次性抛射」对应的源项表达不一样前者是恒定源强后者是源强随时间衰减的脉冲源。我的做法是把所有已知条件列成一张表明确哪些是直接量、哪些需要换算、换算系数出自哪条假设这一步写进论文的假设部分评委基本挑不出毛病。还有一类题目会给出烟囱排气的流速、温度和直径并声明核素活度浓度这时源强可以直接相乘Q 排气流量 × 活度浓度。排气流量等于烟囱截面 × 排气速度注意把工况流量修正到实际温度和压力。4.2 扩散参数表Briggs系数与P-G等级的对应扩散参数是整个模型里唯一能通过查表确定的量。竞赛气象数据往往简化为一行字因此要先把天气描述翻译成 Pasquill 等级。白日按太阳辐射强弱分三个档位夜间按云量分两个档位表里风速行确定最终分类。地面风速(m/s)白天强日照白天中等日照白天弱日照阴天夜间云量4/8夜间云量4/82AA-BBD(F)F2~3A-BBCDEF3~5BB-CCDDE5~6CC-DDDDD6CDDDDD夜间那两列括号代表很稳定的晴夜实际竞赛里出现频率低但一旦出现F 类对应的垂直扩散参数极小地面浓度会非常敏感。我见过有队伍因为取了 E 类而低估了某稳定夜间的浓度和评委参考答案差了近五倍这个细节在论文里要写清楚选取依据。等级确定后用 Briggs 参数表。Briggs 原始公式里 x 单位为 kmσ 单位为 m使用前务必统一量纲。论文里我建议把换算过程写进附录形成一个完整文件清单式的参数说明方便评委查证。一般做法是把上表复制进模型说明并标注「x≤10 km 适用」超出范围时用分段函数处理。4.3 模型检验与灵敏度分析别只给一张图判断模型好坏的指标竞赛中最常用的是归一化均方误差 NMSE 和分数偏差 FB。NMSE 的定义是实测与模拟差的平方除以两者乘积值越小越好FB 反映系统偏差正值表示模拟低估、负值表示高估。没有实测数据时用不同稳定度等级下的模拟结果互相比做敏感性分析。灵敏度分析的标准做法是单因子扰动固定其他参数只把源强、风速、稳定度等级、有效源高分别增减一定比例观察最大地面浓度及其出现距离的变化幅度。通常你会发现风速和源强对浓度的影响最直接而稳定度等级和有效源高决定最大浓度落点的远近这个结论写进论文可以作为「应优先保障气象观测准确性」的论据。模型检验的另一个有效手段是距离衰减曲线。取中心线y0, z0的浓度随下风向距离 x 的变化和理论趋势对比——连续点源在高斯烟羽中近距离浓度随距离近似按 x^-1 衰减远距离受混合层反射影响衰减更快。如果画出来的曲线出现某段异常抬升优先检查扩散参数是否突变或镜像源叠加是否溢出。4.4 一个完整的参数配置清单竞赛团队分工时我习惯给每个人发一张参数清单避免各算各的导致结果对不上。这张表覆盖了从气象到源项的全部输入。参数符号取值示例来源量纲核素释放率Q1.16e10堆芯总活度×份额/时长Bq/s烟囱物理高度hs120题目给定m烟羽抬升高度Δh45Holland公式计算m有效源高H165hsΔhm源高处风速u4.0风廓线幂律换算m/s稳定度等级PGDPasquill查表-衰变常数λ1e-6核素半衰期换算1/s混合层高度zi1000探空或题目假设m这张表的每一行都要写进论文因为参数选取本身就是评分点。特别是烟羽抬升高度很多队伍直接忽略把 H 当作烟囱高度。Holland 公式Δh (vs·d/u)·(1.5 2.68e-3·(Ts-Ta)·d)vs 是排气速度、d 烟囱内径。火电或核电烟囱的热烟气抬升效果显著加与不加 Δh最大落地浓度出现的距离可能从几百米变到几公里这是区分入门队和进阶队的关键细节。5. 避坑指南气体扩散建模里五个容易翻车的细节5.1 现象一地面浓度比预期大一个量级有队伍跑出近地面浓度后总觉得不对劲看起来像是一条浓度很高的狭长带子贴着地延伸。原因是他们把有效源高 H 直接取成了烟囱物理高度忽略了烟气热抬升。原因分析Holland 公式或 Briggs 抬升公式在高温排气场景下会给出几十到上百米的附加高度。H 越大烟羽中心离地面越远地面浓度峰值越低、出现距离越远。忽略抬升等于人为把源压到地面最大地面浓度自然虚高。解决方式烟囱排放参数齐备时先用 Holland 公式计算 Δh再令 H hs Δh。如果题目没给排气温度明确写「假设烟气与环境温度接近不考虑抬升」作为模型假设写进论文不要默默忽略。5.2 现象二同一题目两人跑出完全不同的结果备赛时经常出现 A 同学和 B 同学用同一套公式却得到两倍差异的结果。排查下来十有八九是稳定度等级取挡不一致同一个「晴天三级风」有人取 B 类有人取 C 类。原因分析Pasquill 查表本身就是半经验操作白天中等日照和弱日照的边界很模糊而且 Briggs 参数还有乡村与城市两套系数城市粗糙度更大、扩散更强忘写适用场景会让结果不可比。解决方式全队固定一套判断规则——日照强弱按题目给出的太阳高度角和云量确定并统一使用乡村地形 Briggs 参数。规则写进论文第二章假设同时在参数表里注明稳定度等级的判断依据这样评委复查时结论可复现。5.3 现象三瞬时泄漏场景用烟羽模型浓度永远不衰减如果题目描述的是爆炸性释放某核素瞬间全部进入大气再用连续烟羽的稳态公式去算得出的浓度场只随空间变化而不随时间衰减和物理直觉不符。原因分析烟羽模型的前提是源强在时间上恒定瞬时释放时应采用烟团模型。烟团模型的浓度会随时间增加而摊薄中心浓度随时间按 σy²σz 增长而下降且整体随风漂移。解决方式先把释放过程按时间离散。瞬时释放看作一系列初始时刻不同的独立烟团每个烟团随风输运并扩散空间任一点的浓度是多个烟团贡献的叠加。如果只有一次竞赛保证在论文里明确区分释放类型比追求复杂的多烟团叠加更重要。5.4 现象四等浓度线出现「反向尾巴」统计结果时有人发现等浓度线在下风向某处又向轴上收拢甚至出现与风向相反的拖尾形状看起来像浓度在下风向先增后减又突增。原因分析扩散参数 σz 随距离增长太快当 x 超过 Briggs 公式适用上限通常 10 km垂直扩散参数被高估浓度被稀释得过快。远距离处的模拟浓度偏低曲线在图上呈现非物理的「回升」。解决方式给 σz 设置上限同时把混合层反射项作为约束。当 σz 接近混合层高度说明烟羽已经充分混合垂直方向浓度应趋于均匀而非继续按高斯摊薄。在代码里增加判断若 sz 0.5×h_mix垂直项直接用均匀混合近似这是工程上处理远场扩散的常用做法。5.5 现象五蒙特卡洛不确定性分析不收敛把风速、源强、稳定度作为随机变量做不确定性分析结果却每次运行都不一样置信区间宽得离谱。原因分析随机抽样次数不够或者稳定度等级这种离散变量直接用均匀分布抽样样本点大量集中在物理上不可能的组合上导致结果方差被拉大。更隐蔽的原因是随机种子没有固定两组结果不可复现。解决方式稳定度等级按实际气候频率分布抽样比如中性 D 类占一半以上风速用对数正态分布或威布尔分布源强用对数均匀分布。把随机种子固定样本数至少取 5000 以上。追求效率时用拉丁超立方抽样几百个样本就能覆盖参数空间比纯随机收敛快得多评委也更认可这个方法。6. 从浓度到剂量让模型结果对接应急决策的一个技巧6.1 从浓度场算到吸入有效剂量再画一张等值线图浓度分布画出来只是完成了工程计算评委真正想看的是你能不能回答「哪些区域需要撤离」这个应急决策问题。我的习惯是最后一步做剂量换算把地面浓度场乘上呼吸率和照射时间得到吸入有效剂量再和应急干预限值对比画出剂量等值线。这一步能让模型的输出从「Bq/m³」变成「mSv」语义上直接对接辐射防护论文的价值立刻上一个台阶。剂量换算的公式很简单D C · BR · g · T。C 是地面浓度单位 Bq/m³BR 是呼吸率成年人轻度活动取 1.2 m³/h但竞赛里按连续 24 小时吸入计就需要乘 24g 是吸入单位活度对应的有效剂量转换系数这个系数按核素查标准表竞赛题通常会附在附录里T 是受照时间。换算时要注意单位统一我见过有人把每小时呼吸率和 24 小时代进去忘了乘 24结果低估了剂量。等值线图用 matplotlib 的 contour 画把待评估核素的剂量限值作为分档线画清楚 1 mSv、5 mSv、50 mSv 这几条关键曲线。等值线闭合区域的边界就是撤离范围建议。画完之后我会在图上标出最大剂量点和最远受影响边界并在论文结论处写一段「建议以 XX 等值线作为应急行动分级界线」。这一小段话通常就是评分细则里「模型应用价值」栏的得分点。这个技巧背后其实是一个思维习惯建模的终点不是漂亮曲线而是决策建议。我现在每次带队都会要求先在白板上把「浓度→剂量→限值→撤离距离」这条链写出来再动笔写代码。这样做还有个好处就是能从应急限值反推网格范围和计算精度——如果关心的限值浓度在 20 km 外才出现网格就必须覆盖到 30 km否则等值线在边界处被截断结果很尴尬。希望这个思路在你的竞赛里能派上用场也希望这篇实战拆解能帮到你。本文还有配套的精品资源点击获取