简介一份围绕核电站泄漏后放射性气体浓度分布规律与扩散模型研究的数学建模竞赛完整方案适合参加数学建模比赛的学生、指导教师及从事环境影响评估的研究者参考。方案从连续泄漏与瞬时泄漏两类源出发构建高斯烟羽模型、一维与三维抛物型扩散模型、有限时间泄漏扩散模型以及高斯烟团模型涵盖分离变量法、傅里叶变换、MATLAB求解和福岛核泄漏实际案例分析。内容依据陕西师范大学2011年全国大学生数学建模竞赛模拟赛题整理完整包含赛题原文、模型假设、符号说明、公式推导与结论整理并对无风、有风以及上风/下风不同场景下的浓度预测过程做了细致说明。压缩包共1个doc文档大小1.18MB目录结构清晰便于直接对照学习。已有204人学习下载。1. 核电站泄漏后放射性气体扩散模型数学建模题里最容易被低估的一类近几年的华为杯研究生数学建模竞赛和国赛题目里核安全类问题反复出现核电站泄漏后的放射性气体扩散就是其中典型。这道题表面看是环境工程问题实际上考的是三件事能不能把事故场景抽象成正确的物理模型、能不能把气象参数翻译成模型输入、能不能让评审一眼看出你的浓度分布规律经得起推敲。大多数队伍死在第 2 件事上——高斯烟羽模型公式背得熟但源项怎么定、稳定度怎么查、扩散参数怎么修正全凭感觉。这篇内容把气体扩散模型里最常用的高斯烟羽系做法拆开讲从参数化到数值实现再到反演验证给出一套能直接用到赛题里的技术路径。2. 高斯烟羽模型与核泄漏源项参数化先回答这 3 个问题再套公式2.1 为什么泄漏后的放射性气体扩散首选高斯烟羽系模型放射性气体从烟囱或破口释放后在大气中的输运主要受平均风场支配湍流脉动负责把污染物向周围稀释。对一场持续几小时到几天的泄漏事故污染物在水平方向上的浓度梯度远小于沿风向的变化这种条件下时间平均的浓度场可以用稳态高斯烟羽模型来描述。它的核心假设是污染物的横向和垂向浓度分布都服从正态分布扩散系数随风向距离按幂律增长。这个假设在平坦地形、风速不低于 1 m/s、气象条件相对均匀的工况下误差可控而在数学建模竞赛里题目通常给的就是这种简化场景。相比拉格朗日粒子模型和 CFD 方法高斯烟羽系模型的计算成本低一个量级以上参数都有经验公式可以查且在 10 km 以内的近场浓度预测上精度并不差。对竞赛场景这意味着可以在答辩现场跑参数敏感性、画浓度云图、做源项反演而不是把时间耗在网格收敛性上。常见做法是先建立地基烟羽模型跑通全流程再用烟团模型或 CFD 对局部复杂工况做对比验证这块后文展开。2.2 源项三要素释放率、有效源高与有效清除常数源项参数化是整个气体扩散模型里最影响结果的一步也是竞赛论文里最容易扣分的地方。源项至少包含三个量释放率 QBq/s、有效源高 Hm、有效清除常数 λ1/s。释放率在题目中可能直接给出也可能只给出堆芯破损程度要求自行估算后者需要结合实际泄漏持续时间做归一化把“泄漏总量”换算成“平均释放率”。有效源高不是烟囱的物理高度。高温放射性气体喷出后受浮力作用会继续抬升所以有效源高 H 由物理高度 h_s 加上抬升量 Δh 构成。常见做法是用 Briggs 抬升公式估算中性层结下的抬升量Δh 1.6 · F^(1/3) · u^(-1) · x^(2/3)其中 F 是浮力通量参数m⁴/s³u 是环境风速m/sx 是下风向距离。注意 Briggs 公式在小风速下会出现抬升量异常偏大的情况如果题目给了静风工况需要对该公式做限幅处理。有效清除常数 λ 不是单指放射性衰变常数 λ_r。放射性气体在扩散过程中同时存在干沉积、湿沉积和衰变三者共同构成有效清除率。对惰性气体如 Xe-133、Kr-85沉积项可以忽略λ 近似等于 λ_r对碘和铯同位素沉积项占比显著必须把干沉积速度 v_d 折算进去λ_eff λ_r v_d / H_mix其中 H_mix 是混合层高度mv_d 是干沉积速度m/s取值可以通过题目给的地面粗糙度来定。下表给出竞赛中最常见核素的衰变常数量级方便代入公式做量级判断。核素半衰期λ_r 量级(1/s)常见泄漏场景Xe-1335.2 天约 1.5×10⁻⁶堆芯燃料破损释放Kr-8510.7 年约 2.0×10⁻⁹乏燃料处理I-1318.0 天约 1.0×10⁻⁶安全壳泄漏Cs-13730.2 年约 7.3×10⁻¹⁰燃料严重损毁2.3 从题目文字提取参数的常见坑竞赛题很少把参数表直接摆出来多数情况下藏在“泄漏源高度 xx m”“事故发生在夏季白天”“风速 xx m/s”这类描述里。夏季白天对应不稳定层结夜间多为稳定层结A—F 级的划分直接影响 σ_y、σ_z 的取值。如果题目给了多个时刻的监测浓度数据务必要把释放率 Q 和有效源高 H 当作未知量去反演而不是直接取一个看起来合理的值——后文的第 6 章会专门给出反演方法。此外所有放射性浓度计算都要额外乘一个 exp(-λ_eff · x / u) 的衰变因子这个因子在高风速下可能被低估但在低风速下会明显压制远端浓度影响安全区划定的下风向边界。提取完源项之后下一步就是把风场参数化。3. 扩散参数 σ_y 与 σ_z 的取值稳定度分类、粗糙度修正与混合层顶的封顶处理3.1 Pasquill 稳定度分类与幂律拟合参数表高斯烟羽公式里决定浓度分布形态的核心参数是水平扩散参数 σ_y 和垂直扩散参数 σ_z二者都是下风向距离 x 的函数。工程上最常用的做法是用 Pasquill 稳定度分类法把大气层结分成 A极不稳定到 F稳定共六个等级再套用经验幂律σ_y a · x · (1 b · x)^(-1/2) σ_z c · x · (1 d · x)^(-1/2)其中 a、b、c、d 是经验常数x 单位取 mσ 单位取 m。不同资料手册给出的常数略有差异竞赛中用下面这组常用参数即可评审不会因为常数差 10% 扣分但会因为物理场景和参数不匹配扣分。稳定度ab (1/m)cd (1/m)A0.221.0×10⁻⁴0.200B0.161.0×10⁻⁴0.120C0.111.0×10⁻⁴0.082.0×10⁻⁴D0.081.0×10⁻⁴0.061.5×10⁻³E0.061.0×10⁻⁴0.033.0×10⁻⁴F0.041.0×10⁻⁴0.0163.0×10⁻⁴稳定度等级怎么由题目给的信息确定白天日照强、风速小断定 A 或 B 级阴天且有中等风速时取 D 级夜间晴朗微风取 F 级。一个非常实用的判断是200 m 以上的高架源在 D 级条件下地面浓度通常在下风向 3~6 km 处达到峰值如果题目给的监测点峰值出现在 1 km 以内那稳定度很可能不是 D 级而是 A 到 C 级。这类逻辑在竞赛论文里写一句比大段推导更能体现建模水平。3.2 粗糙度修正与地面反射项的处理幂律参数表是基于开阔平坦地形得到的。题目里如果出现城市、森林、丘陵等粗糙下垫面σ_z 会被压抑修正思路是把 σ_z 乘一个小于 1 的系数同时把σ_y 适当放大。粗略做法是城市地形 σ_z 乘 0.6~0.7σ_y 乘 1.2~1.3这个经验范围在 5 km 内可接受。地面反射项的处理则容易踩坑高架源释放的污染物到达地面后会反射回大气公式里表现为两个高斯项相加exp[-(z - H)² / (2σ_z²)] exp[-(z H)² / (2σ_z²)]当 z 0 即关心地面浓度时这两项相等叠加效果等于在系数上乘 2所以地面浓度可以简写为C(x, y, 0) Q / (π · u · σ_y · σ_z) · exp(-y² / (2σ_y²)) · exp(-H² / (2σ_z²)) · exp(-λ_eff · x / u)注意不能把两个指数项直接相加必须先把 (z-H) 和 (zH) 代入再取 z0否则会得到 double count 的错误结果。3.3 混合层顶对垂直扩散的封顶效应当 σ_z 增长到与混合层高度 H_mix 同一量级时污染物垂直方向受混合层顶限制高斯分布被截断。继续用无界高斯公式会高估地面浓度。经验判据当 σ_z 1.6 · H_mix 时垂直方向近似均匀混合地面浓度改用封顶公式C Q / (√(2π) · u · H_mix · σ_y) · exp(-y² / (2σ_y²)) · exp(-λ_eff · x / u)这个公式里已经隐含了多次反射的贡献H 项消失了。这意味着一旦稳定度很强、混合层很薄有效源高对地面浓度的主导权会被混合层顶取代。竞赛里如果题目给了逆温层高度或风速廓线数据一定要做这个判断它往往是区分优秀论文和普通论文的分水岭。实现层面把 σ_z 1.6·H_mix 作为切换条件写进浓度计算函数比在公式层面硬推更可控。4. 用 Python 在网格上数值化放射性气体浓度场从公式到云图4.1 高斯烟羽地面浓度的最小实现把上一章的公式封装成可复用函数是整道赛题的地基。下面这个 Python 实现基于 numpy输入为网格坐标、源项和气象参数输出为地面浓度网格可直接用于绘图和后续反演。import numpy as np def pasquill_sigma(x, stabilityD): Pasquill-Gifford 幂律扩散参数x 为下风向距离(m)。 params { A: (0.22, 1.0e-4, 0.20, 0.0), B: (0.16, 1.0e-4, 0.12, 0.0), C: (0.11, 1.0e-4, 0.08, 2.0e-4), D: (0.08, 1.0e-4, 0.06, 1.5e-3), E: (0.06, 1.0e-4, 0.03, 3.0e-4), F: (0.04, 1.0e-4, 0.016, 3.0e-4), } a, b, c, d params[stability] # 用带分母修正的幂律形式避免远端发散 sigma_y a * x / np.sqrt(1.0 b * x) sigma_z c * x / np.sqrt(1.0 d * x) return sigma_y, sigma_z def plume_ground(x, y, Q, H, u, lam, stabilityD, H_mixNone): 计算地面(z0)的放射性气体浓度x,y 为基于风向坐标系的坐标。 x np.atleast_1d(np.asarray(x, dtypefloat)) y np.atleast_1d(np.asarray(y, dtypefloat)) sigma_y, sigma_z pasquill_sigma(x, stability) if H_mix is None or np.all(sigma_z 1.6 * H_mix): # 无界或未触顶使用全反射高斯烟羽公式 coeff Q / (np.pi * u * sigma_y * sigma_z) term_y np.exp(-y**2 / (2.0 * sigma_y**2)) term_h np.exp(-(H**2) / (2.0 * sigma_z**2)) decay np.exp(-lam * x / u) else: # 垂直均匀混合混合层封顶公式 coeff Q / (np.sqrt(2.0 * np.pi) * u * H_mix * sigma_y) term_y np.exp(-y**2 / (2.0 * sigma_y**2)) term_h 1.0 decay np.exp(-lam * x / u) return coeff * term_y * term_h * decay这段代码的逻辑分三层。先算扩散参数 σ_y、σ_z再按是否触顶选择高斯全反射公式或混合层均匀公式最后乘上衰变项。要注意 x 和 y 必须是以风向为 x 轴的正交坐标系下的坐标如果原始数据给的是经纬度或平面直角坐标必须先做旋转或映射。参数上Q 的单位是 Bq/su 是 m/slam 是 1/s输出的浓度单位是 Bq/m³。H_mix 设为 None 时不做封顶处理适合低空泄漏或近场计算。4.2 用风向旋转把监测点坐标统一到烟羽坐标实际题目给的是场地直角坐标而风向通常不与坐标轴平行。计算前要把所有坐标点绕泄漏源旋转到以风向为 x 轴的标准坐标系。旋转公式是def rotate_to_wind(x, y, wind_dir): wind_dir 为风向角(度)表示风从哪个方向吹来。 theta np.deg2rad(wind_dir - 90.0) x_r x * np.cos(theta) y * np.sin(theta) y_r -x * np.sin(theta) y * np.cos(theta) return x_r, y_r注意 wind_dir 的定义。气象上通常给“风的来向”北风意味着风向角 0 度、污染物向南输送工程上有时直接给“去向”。拿到题目先确认是哪一个否则整张浓度云图会对称翻转。旋转之后所有 x_r 为负的点位于泄漏源上风向浓度按 0 处理即可不需要代入高斯公式。4.3 浓度分布规律的可视化与网格分辨率选择网格间距选多少直接决定计算量和云图平滑度。常见做法是取 100 m 的网格距覆盖下风向 20 km、侧风向 6 km 的范围网格数约 200×60。这个规模在 numpy 下单次计算毫秒级做参数敏感性扫描毫无压力。绘图时用 contourf 画浓度云图并对浓度取对数或设置阈值截断因为放射性浓度跨越多达 6 个数量级线性色标会让近源区一片白、远端一片黑。典型代码import matplotlib.pyplot as plt x_grid np.arange(-2000, 20000, 100.0) y_grid np.arange(-3000, 3000, 100.0) X, Y np.meshgrid(x_grid, y_grid) C plume_ground(X.ravel(), Y.ravel(), Q3.0e15, H80.0, u4.0, lam1.0e-5, stabilityD).reshape(X.shape) plt.figure(figsize(10, 4)) cf plt.contourf(X / 1000.0, Y / 1000.0, np.log10(C 1e-6), levels18, cmaphot) plt.colorbar(cf, labellog10(C / (Bq/m^3))) plt.xlabel(x / km); plt.ylabel(y / km) plt.title(Ground-level concentration field)这里把浓度加了 1e-6 再取对数是为了避免零值区域出现 -inf。观察云图时会发现两个明显规律一是中心线浓度随下风距离先增后减峰值出现在某个特定位置附近二是横向浓度剖面在每个 x 截面处都保持高斯形状。峰值位置和高度对有效源高 H 和稳定度非常敏感这正是下一章敏感性分析要捕捉的对象。5. 参数敏感性分析、模型对比与竞赛论文的复盘路径5.1 控制变量法扫描 Q、H、u、λ 的影响拿到一个能跑的浓度场之后不要急着堆别的模型先把单因素敏感性分析做扎实。以中心线地面浓度峰值 C_max 和峰值出现距离 x_max 为输出指标控制其他变量逐一扫描结果用表格格式呈现是华为杯数学建模优秀论文里最常见的做法。下面是一次典型扫描的示意结果Q3×10¹⁵ Bq/sH80 mu4 m/sD 级稳定度变量变化范围对 C_max 的影响对 x_max 的影响u 风速2~8 m/s风越大峰值越低峰值位置小幅后移H 有效源高40~160 m高度翻倍峰值约降 1 个量级峰值位置明显后移Q 释放率1×10¹⁵~1×10¹⁶线性正比基本不变λ 清除常数1×10⁻⁶~1×10⁻⁴远端浓度被压缩峰值位置略微前移从表格能提炼出三个可写进论文的结论。第一有效源高对峰值浓度是平方级影响因为公式里 H² 出现在指数项中H 翻倍会让 exp(-H²/2σ_z²) 急剧衰减。第二风速 u 的影响不是简单反比因为风速改变后 Briggs 抬升量也随之变化间接影响 H只在固定 H 的前提下才能看到严格反比关系。第三λ 对近场影响弱、对远场影响强这在划定安全距离边界时非常重要。5.2 高斯烟羽、高斯烟团与拉格朗日粒子模型的边界在哪高斯烟羽模型描述连续稳态释放一旦题目变成“事故发生后 10 分钟内泄漏完”这就是瞬时源或有限时长释放烟羽模型的前提出问题需要切换到高斯烟团模型。烟团模型的思路是把释放过程离散成一系列 puffs每个 puff 以风速输送并同时膨胀地面浓度通过对所有 puffs 求和得到。它的公式复杂度比烟羽高但在竞赛中实现起来并不难只需要把烟羽公式中的 x 替换为每个 puff 中心位置到计算点的距离并对时间维度做数值积分。拉格朗日粒子模型LPM在复杂地形和边界层非均匀场景下更准确但需要气象场驱动数据竞赛题通常不提供所以它的角色是“用来说明模型的局限”而不是实际计算主力。CFD 同理用来做城区绕流和建筑物尾流效应这类局部场景的定性讨论。论文里如果能在灵敏度分析之后补一段“模型适用性对比”说明什么时候烟羽足够、什么时候必须升级到烟团或粒子模型整套建模逻辑就闭环了。5.3 “浓度分布规律”到底要讨论哪些规律竞赛题名里的“浓度分布规律”不是空话它要求从计算结果里归纳出不依赖具体参数的一般性结论。常见做法是画三条曲线一是中心线地面浓度随下风向距离的变化曲线归一化坐标 C·u/Q二是不同稳定度级别下的侧风向浓度剖面验证高斯形状的“胖瘦”差异三是同一稳定度下不同有效源高对应的峰值位置包络线这个往往能和解析解对应上。这三条曲线一旦做出来后半篇论文就有了锚点后续的污染影响评估、安全距离划定都从这个结论展开。6. 用监测点浓度反推泄漏源强一个有资格写进“创新点”的技巧竞赛题通常会给几组监测点浓度数据而 Q 和 H 恰好未知这时把正问题反着用就是反演。做法是从假设的 Q、H 出发算监测点浓度与观测值比对用最小二乘不断修正假设值。这个环节不需要复杂的机器学习scipy.optimize.least_squares 足够用。from scipy.optimize import least_squares # 假设观测点坐标和浓度已知 obs_points [(800, 120), (2000, -250), (3500, 80)] obs_concentration [4500.0, 620.0, 180.0] # 单位 Bq/m3示例数据 def residual(theta): Q_est, H_est theta pred [] for (xi, yi) in obs_points: c plume_ground(np.array([xi]), np.array([yi]), QQ_est, HH_est, u4.0, lam2.0e-6, stabilityD) pred.append(c[0]) return np.log10(np.array(pred)) - np.log10(obs_concentration) res least_squares(residual, x0[1.0e14, 60.0], bounds([1.0e10, 10.0], [1.0e19, 300.0])) Q_opt, H_opt res.x print(fEstimated Q {Q_opt:.3e} Bq/s, H {H_opt:.1f} m)这里对预测值和观测值都取了对数再求残差原因是浓度跨数量级分布如果直接在原始尺度上做最小二乘近源点会主导损失函数而远端信息被稀释。注意残差的构造也把 Q、H 的物理范围写进了 bounds防止优化器给负源强或不合理的大源高。实际使用中监测点数量少于 3 个时反演问题是不适定的常见处理办法是先根据浓度峰值的横向展宽估 σ_y从而反推稳定度级别把稳定度固定后再反演 Q 和 H。验证反演结果有一个盲测技巧把模拟得出的浓度场删掉 20% 的监测点用剩余 80% 反演再去和被删点比对这叫交叉验证。竞赛论文里放一组交叉验证残差图说服力远大于“结果符合预期”这句话。反演失败时优先检查位移观测点坐标是否旋转到风坐标系、观测时刻是否对应稳态释放、监测点在 x 方向上是否落在泄漏源上风向。坐标旋转错了最简单的最小二乘也会给出荒谬的源强这种错误在历届参赛论文里反复出现。反演结果出来后别忘了做不确定性传播。把 Q、H、u 的置信区间代入正问题得到浓度的置信带这个置信带直接决定安全区划定的保守程度。如果评审质疑你的 Q 反演值太大把敏感性表和置信带一起展示比口头解释有说服力得多。本文还有配套的精品资源点击获取