简介一份面向全国大学生数学建模竞赛的核电站泄漏放射性气体扩散建模论文方案重点研究不同泄漏源类型、风速风向等因素下的浓度分布规律。资源包含基于高斯烟羽模型、一维抛物型扩散模型、三维空间扩散模型、有限时间泄漏扩散模型及高斯烟团模型的完整建模推导并结合福岛核泄漏实例分析对我国东海岸和美国西海岸的影响适合数学建模参赛者、环境风险评估人员参考。包体为单个doc文档体积约1.18MB内容涵盖问题重述、模型假设、符号说明、公式推导与MATLAB求解结果可作为竞赛论文撰写与模型复现的底稿。已有204人学习浏览对希望快速掌握气体扩散建模思路或需要完整赛题论文范例的读者具有实用参考价值。1. 核电站泄漏气体扩散模型一份能直接复现的数学建模竞赛全案做数学建模竞赛的人迟早会遇到气体扩散这道坎核电站泄漏、化工厂毒气泄漏、煤矿瓦斯扩散换着花样考同一个物理过程。这份陕西师范大学2011年模拟赛论文恰好把整条技术路线完整走了一遍从连续源与瞬时源的判定到高斯烟羽模型、一维抛物型方程、三维扩散模型、有限时间叠加模型再到有风条件下的模型Ⅴ和高斯烟团模型Ⅵ最后落到福岛核泄漏案例的浓度估算。它覆盖了数学建模竞赛里气体扩散题的完整链条同时有明确的适用边界——瞬时源假设、风向确定假设、大气稳定度等级划分。对正在备赛的队伍来说这篇论文的价值不只是公式而是每个模型的推导路径、参数来源和坑的位置。下文按复现顺序拆开讲。2. 连续源与瞬时源先判定泄漏类型再选模型高斯烟羽模型Ⅰ只在这条路上成立2.1 判定引理三个区间决定三种建模路线论文开篇先做了一个很关键的判定这件事直接决定后面走哪条建模路线。引理给出的是设气体全部气化或泄漏的时间为 T原文未显式定义该符号通常理解为泄漏持续时间平均风速为 u顺风扩散系数为 σ_x当 T·u / σ_x 远大于 1 时泄漏源作为连续源考虑远小于 1 时作为瞬时源考虑介于中间时最难处理、不做讨论。这里没有用严格的量纲分析而是工程上常用的经验判据2011 年竞赛场景下这个处理是够用的。判定区间背后的物理直觉是污染物在空间里铺开的特征时间如果远大于泄漏时间泄漏过程可以压缩成一个瞬间释放的“团”反过来泄漏持续时间很长污染物在空间里已经形成稳定的浓度尾巴就适合用稳态烟羽来描述。实际复现时泄漏时间数据往往拿不到精确值常见做法是先按最不利情况判定为瞬时源因为连续源模型求解更复杂且很多竞赛题里泄漏时间确实远小于扩散时间。2.2 高斯烟羽模型Ⅰ有效源高的处理与地面浓度公式在判定为连续源的前提下论文建立了高斯烟羽模型Ⅰ。核心参数是有效源高 H它由排放口有效高度 h 和抬升高度 Δh 两部分组成抬升高度用 H h Δh 计算。Δh 的来源是热气体的浮升力和初始冲力竞赛里如果没有具体数据一般按排放口高度加上一个经验估值。坐标系的取法是点源在地面的投影为原点z 轴垂直向上有效源位于 (0,0,H)。无风时浓度公式为C(x,y,z) Q / (2π·u·σ_y·σ_z) · exp(-y²/(2σ_y²)) · [exp(-(z-H)²/(2σ_z²)) exp(-(zH)²/(2σ_z²))]这个公式里没有出现 x 的指数项原因是无风稳态条件下 x 方向浓度直接用源强除以风速和扩散参数表示。实际使用中真正被反复调用的是地面浓度公式令 z 0 得到C(x,y,0) Q / (π·u·σ_y·σ_z) · exp(-y²/(2σ_y²)) · exp(-H²/(2σ_z²))注意地面浓度公式里系数是 π 而不是 2π因为地面反射项和直接项叠加后消掉了分母的 2。这个细节很多人抄公式时会抄错导致浓度值直接翻倍。令 y 0 得到沿 x 轴的浓度分布。进一步对 x 求导并令导数为零可以求出地面浓度最大值出现的位置但 σ_y、σ_z 都是 x 的函数手算很麻烦常见的做法是数值扫值。2.3 有风时与无风时的差别以及参数 σ 的取值依据有风且风向稳定时高斯烟羽模型的形式变为C(x,y,z) Q / (2π·u·σ_y·σ_z) · exp(-y²/(2σ_y²)) · [exp(-(z-H)²/(2σ_z²)) exp(-(zH)²/(2σ_z²))]与无风时比对差异在于分母中出现了风速 u且 x 方向由坐标位置直接决定。σ_y、σ_z 的取值是关键论文采用的是 Pasquill-Gifford 扩散曲线法这种方法查图得到的是下风向距离 x 对应的扩散参数工程上更常用的是幂函数形式σ_y a·x^b σ_z c·x^d其中系数 a、b、c、d 按大气稳定度等级和地面粗糙度查表获得。要注意 P-G 曲线适用于平坦均匀地形、扩散距离在几百米到几十公里的范围超出这个范围外推误差会显著放大。论文里选取的平均风速为离地面 10 米高度处的值这也是气象部门的标准观测高度复现时要保持同一口径。大气稳定度等级是另一个容易出错的地方。论文表 1 给出六级别划分A 极不稳定、B 不稳定、C 弱不稳定、D 中性、E 弱稳定、F 稳定划分依据是地面风速和太阳辐射强度强、中、弱以及云量。白天强辐射加小风速对应 A 级夜间多云加较大风速对应 D 级或 E 级。风速越大大气越接近中性 D 级这个直觉可以用来检验查表结果是否合理。提示P-G 扩散曲线法的原始数据是图直接从图上读数精度有限。复现竞赛题时建议用幂函数拟合值并把拟合来源写进附录阅卷时更好交代。3. 无风瞬时泄漏三件套一维抛物型模型Ⅱ、三维模型Ⅲ、有限时间叠加模型Ⅳ3.1 一维抛物型扩散模型Ⅱ从微元分析到分离变量法判定为瞬时源之后论文先建立最简单的一维模型。假设气体只沿 x 轴直线扩散取核电站为原点设 C(x,t) 为 t 时刻 x 处的浓度。推导用的是微元分析法考虑 x 到 xΔx 段Δt 时间内浓度变化量等于流入量减流出量。流入流出量由傅里叶实验定律给出即单位时间通过单位面积的流量与浓度梯度成正比方向从高浓度向低浓度比例系数就是扩散系数 k。整理后得到一维抛物型方程∂C/∂t k·∂²C/∂x²求解用分离变量法令 C(x,t) X(x)·T(t)代入后拆成两个常微分方程。时间部分解的形式是 T(t) ∝ exp(-λ²kt)空间部分是 X(x) ∝ A·cos(λx) B·sin(λx)λ 无边界条件限制可取任意实数所以最终解是这些基解的叠加对应积分形式的通解。论文强调了 λ 0 时解会随时间无限增长物理上不合理因此舍去只保留衰减模式。这里要提醒一个坑λ 的正负符号处理不当会直接让浓度发散。分离变量时时间方程取 -λ²k 而不是 λ²k这点上课时候老师讲得轻描淡写自己做一遍很容易卡住。一维模型的适用场景很窄实际放射性气体必然往三维空间扩散它更多是作为教学铺垫存在真正的价值在于展示扩散方程的标准推导流程。3.2 三维抛物型扩散模型Ⅲ高斯公式与质量守恒的联立三维模型Ⅲ是整篇论文的枢纽。建立以核电站为原点的三维坐标系z 轴铅直向上浓度记为 C(x,y,z,t)。推导思路仍是质量守恒取任一封闭曲面 S围成区域 ΩΔt 时间从 S 外进入 Ω 的放射性气体质量等于 Ω 内浓度变化引起的质量增量。按扩散定律单位时间通过单位法向面积的质量通量是 -D·∇C其中 D 是扩散系数张量主对角元分别为 D_x、D_y、D_z。对封闭曲面做通量积分利用高斯公式把面积分换成体积分再和浓度变化项联立就得到三维扩散方程∂C/∂t D_x·∂²C/∂x² D_y·∂²C/∂y² D_z·∂²C/∂z²初始条件取 C(x,y,z,0) Q·δ(x)·δ(y)·δ(z)即瞬时点源边界条件是地面不吸收放射性气体对应 z0 处的反射边界。这个边界条件会体现在解的镜像项上。求解用傅里叶变换法把三维方程分解成三个一维方程先对空间变量做傅里叶变换解出时间演化再做逆变换得到C(x,y,z,t) Q / (8·(π·t)^(3/2)·(D_x·D_y·D_z)^(1/2))·exp(-(x²/(4D_x·t) y²/(4D_y·t) z²/(4D_z·t)))这个解的形状是三维高斯分布等浓度面是椭球面椭球的扁率由 D_x、D_y、D_z 的比值决定。扩散系数各向同性时退化为球对称高斯分布。论文中地面不吸收气体的条件在解里体现为 z 用 |z| 处理或在 z0 处加倍复现时建议直接按地面反射处理否则地面浓度会低估一半。3.3 有限时间连续排放模型Ⅳ瞬时源叠加的思路实际问题里核泄漏不是瞬间放完而是持续一段时间。论文的处理方式很朴素把 [0, T] 时间内的连续排放看成无数个瞬时源的叠加每个瞬时源在 τ 时刻释放的量为 Q(τ)到 t 时刻对空间点 (x,y,z) 的浓度贡献为三维高斯分布把所有 τ 从 0 到 T 积分C(x,y,z,t) ∫₀^T Q(τ) / (8·(π·(t-τ))^(3/2)·(D_x·D_y·D_z)^(1/2))·exp(-(x²/(4D_x·(t-τ)) y²/(4D_y·(t-τ)) z²/(4D_z·(t-τ))))·dτ源强恒定 Q(τ) Q 时这个积分可以用误差函数 erf 表达论文给的是含 erf 的闭式结果。误差函数的定义是 erf(x) (2/√π)·∫₀ˣ exp(-s²) dsMATLAB 里直接用 erf 函数调用不需要自己算查表值。这个模型的物理意义泄漏停止后浓度场不是立刻消失而是按高斯分布继续扩散摊薄。t 远大于 T 时模型Ⅳ的结果趋近瞬时源模型Ⅲ这也反过来验证了瞬时源判定区间的合理性。% 模型Ⅳ的浓度计算示例恒定源强T 秒内连续排放 % 输入参数源强 Q0扩散系数 Dx Dy Dz泄漏时长 T观测时刻 t观测点坐标 x y z function C model4_continuous(Q0, Dx, Dy, Dz, T, t, x, y, z) if t 0 C 0; return; end tau_upper min(T, t); % 积分上限不能超过 t也不能超过泄漏时长 if tau_upper 0 C 0; return; end % 数值积分每个瞬时源的贡献叠加 tau linspace(0, tau_upper, 2000); dtau tau(2) - tau(1); Csum 0; for i 1:length(tau) dt t - tau(i); if dt 0, continue; end expo -(x^2/(4*Dx*dt) y^2/(4*Dy*dt) z^2/(4*Dz*dt)); Csum Csum exp(expo) / ((pi*dt)^1.5 * sqrt(Dx*Dy*Dz)) * dtau; end C Q0 * Csum; end这段代码里有两个必须注意的点。积分上限取min(T, t)泄漏还没结束的时刻 t T只有前 t 秒的源有贡献泄漏已结束的时刻 t T积到 T 即可。高斯核里的dt t - tau如果为零指数项分母为零会爆炸所以做了防御性跳过。用 2000 个网格点做数值积分对竞赛场景精度够用如果要把这个代码交到论文附录建议改用integral函数配合匿名函数做自适应积分。注意代码里 D_x、D_y、D_z 的量纲必须是 m²/s。竞赛题给的数据往往混着 km²/h 之类的单位换算错误会让浓度结果差好几个数量级这是最隐蔽的坑。4. 有风条件下的扩散模型Ⅴ与高斯烟团模型Ⅵ上风下风怎么算、参数怎么定4.1 模型Ⅴ的推导微小长方体流量平衡有风且风向确定时论文建立了模型Ⅴ。坐标系取风向为 x 轴正方向z 轴垂直向上原点在核电站地面投影。推导核心是对微小长方体做流量平衡x 方向流量包含两部分一是风速 u 带来的平流项 u·C二是分子扩散项 -D·∂C/∂xy、z 方向没有平流只有扩散项。单位时间内 x 方向流入减流出的净量是 -∂(u·C)/∂x D_x·∂²C/∂x²考虑风速恒定且不可压缩平流项简化为 -u·∂C/∂x。三个方向合并后得到对流扩散方程∂C/∂t u·∂C/∂x D_x·∂²C/∂x² D_y·∂²C/∂y² D_z·∂²C/∂z²这里的 u 是环境风速在 x 方向的分量风向与 x 轴正方向一致时取正相反时取负。论文推导过程中强调了截面积和流入流出的对应关系这是标准的对流扩散方程建模路径工程上大量使用。求解的思路是坐标变换。令 x x - u·t即把坐标系跟着烟气团一起移动对流项就消掉了方程退化回无风扩散方程的形式。初始条件仍为瞬时点源解出来之后把 x 换回 x - u·t得到带平移的三维高斯分布。t 时刻气云团中心在 (u·t, 0, 0)浓度分布为C(x,y,z,t) Q / (8·(π·t)^(3/2)·(D_x·D_y·D_z)^(1/2))·exp(-((x-u·t)²/(4D_x·t) y²/(4D_y·t) z²/(4D_z·t)))地面浓度取 z0 时注意反射项处理。论文给出了地面浓度公式⑺复现时如果只用 z0 代进去而忽略地面反射下风处的浓度会偏低因为实际地面反射会让近地面浓度比自由空间同位置高一倍。4.2 高斯烟团模型Ⅵ正态分布假设与 P-G 扩散参数模型Ⅵ是另一种思路直接假设浓度场在空间上服从正态分布不需要从扩散方程推起。公式形式是三个方向各自带扩散参数 σ_x、σ_y、σ_z 的高斯分布乘积源强 Q 放在系数里C(x,y,z,t) Q / ((2π)^(3/2)·σ_x·σ_y·σ_z)·exp(-((x-u·t)²/(2σ_x²) y²/(2σ_y²) z²/(2σ_z²)))扩散参数 σ 不再是常数而是随下风向距离和时间增长的函数。论文采用 Pasquill-Gifford 扩散曲线法确定 σ 值方法本身来自大量野外示踪实验适用于平坦地形、近地面释放、扩散距离不超过几十公里的场景。σ 与大气稳定度的关系是核心稳定度从 A 到 Fσ_y、σ_z 的增长率依次降低F 级最稳定时扩散最慢。水平扩散参数和垂直扩散参数在 P-G 曲线里是两组不同的图水平方向用 σ_y、垂直方向用 σ_z。实际工程中σ_z 在稳定条件下有一个明显的上限因为大气边界层顶抑制垂直扩散这个效应 P-G 曲线有体现但超出边界层高度后模型失效。模型Ⅴ和模型Ⅵ在物理上是等价的高斯烟团模型的解其实就是扩散系数恒定条件下对流扩散方程的解析解只是扩散参数的表达方式不同。论文里同时建立两个模型目的就是互相验证。复现时建议两个模型都用然后在论文里加一个相对误差对比表这种做法竞赛评阅时比较容易拿分。4.3 上风处与下风处风速符号变换的具体操作问题三要求预测上风和下风 2 公里处浓度。下风处直接套模型公式风速取正值上风处由于污染物是从源点逆风方向扩散过来的对应的是风向坐标里 x 为负的位置。实际操作上不用重推公式把下风处公式里的 x 取负数、u 取正值或者把 u 取负值、x 保持正数数学上等价。论文明确说“只需将风速变为 -u 即可”说的就是这个变换。但有一个物理细节容易忽略上风处能测到浓度靠的是分子扩散和湍流扩散它们的作用方向与风向相反强度远弱于平流。所以上风处的浓度通常比同距离下风处低几个数量级。如果算出来的上风浓度和下风浓度差别不大说明参数设置可能有问题常见的情况是扩散系数设得过大或者风速设得过小。% 模型Ⅴ和模型Ⅵ的地面浓度对比计算 % 参数源强 Q风速 u扩散系数或扩散参数观测点坐标 x,y Q 1e6; % 源强 Bq/s具体量级按题目给的数据调整 u 5; % 风速 m/s离地 10m 观测值 Dx 50; Dy 40; Dz 20; % 扩散系数 m^2/s % 下风 2km 处模型Ⅴ地面浓度 t 2000 / u; % 气云团到达下风 2km 的时间 x_down 2000; y_down 0; z_ground 0; C5_down Q / (8*(pi*t)^1.5*sqrt(Dx*Dy*Dz)) * ... exp(-((x_down-u*t)^2/(4*Dx*t) y_down^2/(4*Dy*t) z_ground^2/(4*Dz*t))); % 上风 2km 处x 取负风速符号处理 x_up -2000; C5_up Q / (8*(pi*t)^1.5*sqrt(Dx*Dy*Dz)) * ... exp(-((x_up-u*t)^2/(4*Dx*t) y_down^2/(4*Dy*t) z_ground^2/(4*Dz*t))); fprintf(模型Ⅴ下风 2km 浓度 %.6g Bq/m^3\n, C5_down); fprintf(模型Ⅴ上风 2km 浓度 %.6g Bq/m^3\n, C5_up);上风和下风浓度都用同一个到达时间 t 2000/u这是刻意简化的做法因为上风处浓度到达峰值的时间与下风处并不相同严格计算分别取各自的到达时间。竞赛中两组结果做对比分析时时间坐标对齐会让曲线更好看但如果要谈物理意义这个处理不够严谨需要在论文里说明。大气稳定度典型场景σ_y 增长趋势P-G 曲线适用性A极不稳定强日照 小风速快适用但外推谨慎B不稳定中等日照 小风速较快适用C弱不稳定中等日照 中风速中等适用D中性阴天或大风较慢最常用数据最全E弱稳定薄云夜晚慢低矮源适用F稳定晴朗夜晚最慢慎用垂直扩散受限5. 避坑记录复现这套核泄漏气体扩散模型时最容易翻车的地方5.1 现象浓度算出来比参考值大好几倍检查公式才发现系数不对原因无风高斯烟羽模型的地面浓度公式反射项叠加后系数从 2π 变成 π。网上流传的版本经常抄漏这个变换直接把有风公式里的 z-H 和 zH 两项相加z0 时得到系数 2/2π但正确的地面浓度公式是 1/π。同样的错误也会出现在建模软件里有人把模型Ⅴ、Ⅵ写成了系数固定版本。解决写代码前先把公式推导一遍确认反射项合并后的系数变化。批量计算时统一用一个函数封装浓度公式不要散落在各个脚本里。校对方法是把 z 设得很大此时反射项影响可以忽略两个版本算出来的值应当一致。5.2 现象扩散系数的数量级明显偏离物理常识比如 D 值取到 10⁵原因题目给的气体扩散速度和泄漏速度往往不是标准单位制直接把 km/h 当成 m/s 代入公式或者把扩散系数与风速混为一谈。扩散系数 D 的量纲是 m²/s典型大气扩散条件下数量级在 10⁰ 到 10² 之间10⁵ 这个量级意味着瞬时扩散到全城明显不合理。解决所有输入数据先统一换算到 SI 单位制再进公式。建立一个参数检查表把每个符号的量纲标注清楚代入前自查一遍。用 MATLAB 算完后用数量级粗检——放射性气体浓度结果如果在 10⁻³ 以下或 10⁶ 以上优先怀疑单位换算而不是模型本身。5.3 现象高斯烟团模型算出来上风浓度和下风浓度差不多曲线几乎重合原因扩散系数设置过大平流项被淹没或者风速 u 设置过小污染物在采样时间内已经逆风扩散到上风很远。实际大气扩散中平流占主导上风浓度应该显著低于下风。解决先确认风速取值是否在 P-G 曲线适用的 1-10 m/s 范围内再检查扩散系数是否超过正常大气量级。如果数据本身没问题把上风和下风的到达时间分别计算不要共用一个 t。正常结果下上风 2 km 的浓度应该比下风同距离低一到两个数量级。5.4 现象查 P-G 曲线时取错了大气稳定度等级E、F 级用在下风远距离浓度分布明显偏窄原因P-G 曲线是按实验数据外推的F 级条件下 σ_z 增长缓慢需要地面粗糙度和边界层高度配合。很多人在远距离外推时不看 σ_z 上限算出的浓度在几十公里外还保持窄高分布这与实际观测不符。解决使用幂函数拟合参数时限定适用距离范围。论文的福岛案例里扩散距离超过数百公里已经超出 P-G 曲线的适用范围这时需要说明这是粗估结果只能看数量级不能抠精确值。在论文里写上“P-G 曲线外推受限”比假装精确要诚实得多评阅时也更好交代。5.5 现象福岛案例里初始浓度和风速都是模拟值结果画出来曲线形状没问题但绝对数值对不上新闻里的监测数据原因真实核事故的源项参数源强、释放持续时间、有效源高几乎不可能精确获取模拟值只能是量级猜测。论文里也承认“具体数值是模拟给出的与实际情况存在偏差”。这不是模型错是输入数据错。解决把源项标成“模拟参数”并在灵敏度分析里给出浓度随风速变化的表。这既是论文的加分项也体现了对模型边界条件的认知。灵敏度分析的具体做法下一章展开。6. 福岛案例的还原技巧参数模拟、风向改变的坐标变换和灵敏度分析三招验证模型6.1 浓度数量级怎么对上先定到达时间再反推源强福岛案例的参数设定是典型的逆向思维。论文已知日方到我国东海岸约 2000 km风速取 5 m/s到达时间约为 4.6 天。先用到达时间把扩散参数 σ_y、σ_z 从 P-G 曲线查出来再根据监测到的浓度数量级反推源强 Q。这个方法在真实事故应急里也常用——源强未知时用远处监测浓度反算源项。反推源强时要注意模型对源强是线性响应浓度正比于 Q所以用一个点测值就能定标。选定参考点后其他点的浓度预测完全由模型决定不再有自由参数。这是验证模型的最佳方式用上海监测点反推 Q用洛杉矶的浓度做预测对比。6.2 风向改变的处理坐标旋转矩阵是唯一要补的数学工具模型的局限在于假设风向全程不变。实际福岛核泄漏期间风向多次变化论文在模型改进部分给出了单次风向改变的坐标变换方法。思路是建立基准坐标系和风向坐标系风向改变前后做坐标旋转。基准坐标系为 (X,Y)风向坐标系为 (x,y)两坐标系夹角为 θ 时坐标关系为x X·cosθ Y·sinθ y -X·sinθ Y·cosθ或者反过来 X x·cosθ - y·sinθ取决于旋转方向约定。将变换关系代入地面浓度公式就能得到基准坐标系下任意时刻任意点的浓度。风向改变多次时用分段函数在时间轴上拼接每个时间段用各自的风向参数。实际操作中风向改变次数多于两三次后解析表达式会变得非常臃肿。常见做法是转成数值模拟把时间步长取小每个时间步按当前风向坐标变换计算浓度场的平移和扩散。MATLAB 里用矩阵运算做这个变换效率很高几百个时间步几秒钟算完。6.3 灵敏度分析作为模型验证和写作技巧论文末尾的灵敏度分析是一个很重要的写作技巧固定初始浓度和其他参数让风速从 1 m/s 变到 3.5 m/s浓度从 1.2428 降到 0.9228灵敏度从 4.9% 上升到 25.7%。表格里每一列算了一次浓度还附了灵敏度公式。这个分析直观地告诉读者风速是模型里最敏感的参数之一风速越大、稀释越快浓度越低。更关键的是它让论文从“求解一个问题”上升到“理解一个问题”的层面竞赛评阅时这是加分项。复现时把灵敏度表里的数据自己跑一遍效率和正确性都能得到验证。从那以后我每次复现这类扩散模型都强制走一遍这三步先做单位制检查、再跑一个参数扫描看灵敏度、最后对比两个独立模型的浓度曲线是否同一个数量级。三个检查都过了结果才敢写进报告。福岛那组模拟数据今天看依然粗糙但整套建模思路——从瞬时源判定到风向变换再到灵敏度分析——放在任何一届数学建模竞赛里都算完整。希望帮到你。本文还有配套的精品资源点击获取