基于高斯模型的放射性物质扩散模型讲解

1、放射性气体扩散浓度预估模型 【摘要】本文是以日本地震引起的福岛核电站的核泄漏为背景， 并以给出的数据 为基础，研究某一假设核电站的核泄漏问题。 我们通过收集相关的资料， 并结合 题目给出的数据， 建立了高斯模型、 连续点源高斯扩散模型解决了题目提出的四 个问题。 针对问题一：考虑到泄漏源是连续、均匀和稳定的，我们运用散度、梯度、 流量等数学概念， 通过“泄漏放射性物质质量守恒”、 “气体泄漏连续性定理”、 Guass公式及积分中值定理得到了无界区域的抛物线型偏微分方程，然后再通过 电源函数解出空间任意一点的放射性物质浓度的表达式， 把此表达式定为模型一 的前身。鉴于放射性物质的扩散受到诸多因素

2、的影响，如：泄漏源的实际高度、 地面反射等。我们以泄漏口为坐标原点建立三维坐标系，通过“像源法”处理地 面反射对放射性物质浓度的影响， 并由此对模型一的前身进行修正完善， 得到模 型一：高斯模型， 即放射性物质浓度的预测模型。 最后我们模拟了放射性物质无 风扩散仿真图。 针对问题二：当风速为 k m/s 时，我们根据放射性核素云团在大气中迁移和 扩散的数值计算的基本方法和步骤，并以泄漏点源在地面的投影点为坐标原点， 以风向方向为 x轴，铅直方向为 z轴，与 x轴水平面垂直方向为 y 轴建立三维坐 标系，地面的反射作用同样利用“像源法”进行处理，得到连续点源高斯扩散模 型。考虑到地面反射、 烟云

3、抬升、 放射性物质自身的沉降及雨水的吸附等对浓度 的影响，我们对连续点源高斯扩散模型进行了修正， 建立了修正的连续点源高斯 扩散模型。最后利用大气稳定度确定了扩散参数，进而求解了模型。 针对问题三：经分析，问题三的提出是以问题二为基础的，模型三的建立只 需要将模型二加以调整即可。我们以风速方向为 x 轴正方向，将风速与放射性物 质的扩散速度进行矢量运算，此问题则转化为求 (L,0, z) 和( L,0, z)两点处的放 射性物质浓度， 由此建立模型三， 即上风和下风 L 公里处放射性物质浓度浓度的 预测模型。 针对问题四： 首先，我们通过网络收集了相关数据， 然后，我们结合模型二、 模型三对数

4、据进行整理代入， 算出了日本福岛核电站泄漏的放射性物质扩散到中 国东海岸和美国西海岸的浓度分别为 4.2429 10 3g/ m3、 2.3854 10 4g/m3。 关键词： 高斯模型 连续点源高斯扩散模型 核泄漏 一 问题的重述 1.1 问题背景 目前，核电站的发展能带来巨大的经济效益和社会效益，但核电站一旦发生 核泄漏，将会给人们的生命健康和周边环境带来巨大的危害性影响。 2011 年 3 月日本的福岛核电站的放射性气体的核泄漏事件更让我们关注放射性气体泄漏 时的浓度问题。 因此，正确的测出大气中放射性物质的浓度在环境监测和安全评 估中具有重要意义 1.2 问题提出 有一座核电站遇自然灾

5、害发生泄漏，浓度为 p0 的放射性气体以速度 m kg/s 匀速排出，在无风的情况下，以速度 s m/s 匀速在大气中向四周扩散 . 1) 在无风的情况下， 建立一个描述核电站周边不同距离地区、 不同时段放射性 物质浓度的预测模型。 2) 当风速为 k m/s 时，给出核电站周边放射性物质浓度的变化情况。 3) 当风速为 k m/s时，分别给出上风和下风 L 公里处，放射性物质浓度的预测 模型。 二 符号说明 i(i x,y,z) 空间任意一点的放射性物质的扩散系数 C(x,y,z,t) 空间任意一点的放射性物质浓度 s 放射性气体的扩散速度 , m/ s Q0 泄漏源泄漏的放射性物质总量 空

6、间域 V 空间域其体积 S 一规则的球面面积 Q1 在 (t,t t) 内通过 的流量 Q2 内放射性物质的增量 t 任意扩散时刻 H 泄漏源距地面的实际高度 H 烟云抬升高度 h 泄漏源有效高度 u 泄漏源高度处的平均风速， m/ s Q 源强， kg / s x, y, z 用浓度标准差表示的 x,y,z 轴上的扩散参数 vs 放射性气体出口流速 ,m/s D 泄漏源出口的有效内径 QH 泄漏源的热排放率， kw Ts 泄漏源出口处温度， K Ta 环境大气平均温度， K Vs 沉降速度， m/ s 空气的动力粘性系数 雨水吸附系数 h0 太阳高度角 三 模型假设 1、扩散过程中浓度在 y

7、、 z轴上的变化分布是高斯分布。 2、放射性物质的扩散看作是空间某一连续点源向四周等强度地瞬时释放放射性 物质，放射性物质在无穷空间扩散过程中不发生性质变化，且不计地形影响。 3、放射性物质扩散服从扩散定律，即单位时间通过单位法向面积的流量与它的 浓度梯度成正比。 4、放射性物质在穿过降雨区域时，其强度由于雨水的吸收而减少，减少比率为 常数。 5、假设地面对放射性气体起全反射作用，地面和海面对放射性气体没有吸附， 将海面视为平原地区 6、假设风向为水平风向，且风向风速不随时间变化。 7、扩散过程中不考虑泄漏点内部温度的变化对气体扩散的影响 四 问题的分析 4.1 问题（ 1）的分析 一座核电站

8、遇自然灾害发生核泄漏， 浓度为 p0 放射性气体以速度 m kg/s 匀速 排除，这近似于放射性物质源是连续均匀稳定的。 在无风情况下， 放射性气体以 速度 s m/s匀速在大气中向四周扩散，放射性气体的扩散服从扩散定律，即单位 时间通过单位法向面积的流量与它的浓度梯度成正比。 在这些条件下， 我们明确 了要研究的问题是 点源连续泄漏的扩散问题，为了使建立的模型更加贴近实际， 需考虑地面反射、核泄漏源的实际高度、降雨等因素对浓度分布的影响。由“扩 散定律”“放射性物质质量守恒定律” “气体泄漏连续性定理” 可得出无界区域的 抛物线型偏微分方程。 再通过假设条件建立未考虑地面反射、 核泄漏源的实

9、际高 度、降雨等浓度影响因素的初步模型， 然后从这些影响因素对模型进行完善， 最 终得出 核电站周边不同距离地区、不同时段放射性物质浓度的预测模型。 4.2 问题（ 2）的分析 本问是探究风速为 km/ s时，核电站周边放射性物质浓度的变化情况。当环 境中空气流动产生风力时， 在均匀湍流场中， 扩散参数与下风向距离的关系是明 确的，核泄漏时间较长时， 可认为扩散是稳定的。 在下风向的湍流扩散相对于风 力引起的移流相可忽略不计， 在流动方向建立 x 轴，不考虑横向速度和垂直速度。 根据假设， 空间中放射性物质的浓度服从高斯分布， 可利用连续点源放射性物质 的高斯扩散模型。 放射性物质在大气中扩散

10、受诸多因素影响， 考虑泄漏源有效高 度、放射性物质自身重力产生的重力沉降、 雨水的吸附等因素对放射性物质浓度 的影响是必要的， 通过这些影响因素对高斯模型进行修正， 然后利用修正后的高 斯模型探究核电站周边放射性物质浓度的变化情况。 4.3 问题（ 3）的分析 本问是要求当风速为 km/ s时，建立上风和下风 L 公里处的放射性物质浓度 的预测模型。经分析，此问是问题（ 2）的延伸，我们只需建立合适的坐标系， 将此问题转化为求具体两处的放射性物质的浓度， 便能得出上风和下风 L 公里处 的放射性物质浓度预测模型。 五 模型的建立与求解 5.1 模型一的建立与求解 5.1.1 模型一的初步建立

11、以核泄漏点正下方的地面为坐标原点 （0,0,0） ，平均风向为 X 轴、指向下风 方向，铅直方向为 Z 轴，水平垂直于风向轴（ X 轴）为 Y 向，建立空间坐标系， 则核电站泄漏点 O距有效地面的高度为 H ，则泄漏点位置坐标为 O（0,0, H ）。 图一：空间坐标系示意图 将气体从泄漏源泄漏时刻记作 t=0，时刻 t 无穷空间中任意一点坐标为 ( x,y,z) 的浓度记为 C(x,y,z,t), 根据假设 2，单位时间通过单位法向面积的流量与浓度 梯度成正比，则： q i gradC( 1) i (i x, y, z) 是扩散系数， grad 表示梯度，负号表示由浓度高向浓度低的地 方扩散

12、。 考察空间域 ，其体积为 V ，包围 的曲面为 S， S为一规则的球面， S外 法线向量为 n (- x,- y ,1) 。则在 (t,t t)内通过 的流量为： zz tt 2) Q1 t q nd dt s 内放射性物质的增量为： 3) Q2C(x,y,z,t t) C(x, y, z, t) dV V 从泄漏源泄漏的放射性物质的总量为： Q0 tt tp0dVdt 根据“质量守恒定律”和“气体泄漏连续性原理” ，单位时间内通过所选曲 面 S的向外扩散的放射性物质的量与 S 曲面内放射性物质增量之和， 等于泄漏源 在单位时间内向外泄漏的放射性物质。则： Q0 Q1 Q2( 5) 即： C

13、(x,y,z,t t) C(x,y, z, t ) dV t q nd dttt tp0dVdt (6) Vs （其中 div 是散度记号 ） 7） 又根据曲面积分的 Gauss公式： q nd divqdV sV C(x,y,z,tt) C(x, y,z,t) tdV divqdVdttt tp0dVdt V t t V tt C C(x,y,z,t t) C( x, y, z, t)t kdiv ( gradC )dt lim lim t t t 0 t t 0 t 由以上两式得： dV t divqdV tp0dV t V t V 即为： Ct dVdivqdVp0dV（ 8） V t

14、V V 由以上公式并利用积分中值定理得： 2 2 2 C C C C i div（ gradC ） i（ 2 y 2 z 2 ）,t 0, x,y,z （9） t x y z 这是无界区域的抛物线型偏微分方程，根据假设1，初始条件为作用在坐标 原点的电源函数，记作 C（ x, y, z, 0） 0Q （x ,y ,z ）（10） Q0表示泄漏源漏泄释放的放射性物质总量， （x, y, z）是单位强度的电源函数。 方程（ 9）满足方程（ 10）的解为： x2 y2 z2 C(x,y,z,t) (4 i t)3/2 4 it ,i (x,y,z) 此模型只是在不考虑风速的情况下建立的，但为了使模型

15、具有更加的实用 性，下面我们将考虑泄漏源的实际高度、地面反射、降雨等因素对浓度的影响， 完善模型。 5.1.2 模型一的最终建立 1. 地面反射对模型的完善 泄漏源有一定的高度，且泄漏点源是连续点源，则泄漏点源可视为高架连续 点源，考虑到地面对扩散来的放射性气体有反射作用，根据假设4，地面对到达 地面的扩散气体完全反射。 这儿可认为地面就像镜子一样， 对放射性气体起全放 射作用，可用“像源法”处理，如图 3，建立三个坐标系，一是以泄漏源（实源） 为坐标原点； 二是以泄漏源在地面的投影点为原点， p点是空间的任意一点， 坐 标为 （ x, y, z） ；三是以泄漏源关于地面的像对称源（像源）为原

16、点。把p点放射 性气体浓度看成两部分（实源与像源）作用之和。 图二： 高架连续点源扩散示意图 从以上分析知， p 点放射性气体的浓度为实源和像源的放射性气体扩散至此 （ z H ） 2 点浓度的叠加。 则实际泄漏源（实源）对 p点的浓度贡献部分可用 e 2 z2 来表示； 因为地面对扩散物质完全反射， 则像对称源（像源）对 p点的浓度贡献部分可用 （z H ） e 2 z2 来表示。于是对（ 11）式所修正完善的模型为： C(x, y, z,t) mp0 (4 i)3/2 t1/2 (z H )2 (z H)2 4 zt z ,i (x, y, z) 12) 我们自己模拟一组数据， 利用 ma

17、tlab进行仿真模拟， 可实现该模型的模拟图 像。假设扩散系数=0.00001，放射性物质的初始浓度C =100，扩散时间 i t=1000000，放射源总量 Q =1000000。（程序见附录一） 图三：核泄漏无风扩散 5.2 模型二的建立与求解 5.2.1 模型二的建立 放射性核素云团在大气中迁移和扩散的数值计算基本上可分为二步。 第一步 根据大气动力学理论进行所关心区域中风场的计算，其理论基础是大气运动方 程、连续性方程、状态方程、热力学方程和水汽方程构成的基本方程组。在大气 科学研究领域中， 已有多个实用的大气环流模式。 第二步进行已知风场中放射性 核素云团迁移和扩散的计算， 可采用类

18、似于处理大气污染的方法， 假设放射性核 素云团不影响大气流体速度和温度，求解放射性核素的连续性方程。 当风速为 k m/ s时，利用连续点源高斯扩散模型分析核电站周边放射性物质 浓度的变化情况。此泄漏点源是有边界点源，泄漏点源的实际高度为 H 。以泄 漏点源在地面的投影点为坐标原点， 以风向方向为 x轴，铅直方向为 z轴，与 x轴 水平面垂直方向为 y 轴建立三维坐标系，由于扩散过程中浓度在 y 、 z 轴上的变 化分布符合高斯分布，所以下风向的任意一点 C(x,y,z) 的浓度函数为： 22 C(x, y, z) A(x)e ay e bz(13) 根据概率统计我们可以得出方差的表达式为：

19、2 y2C(x, y, z)dy 0 0 C(x, y,z)dy 2 0 z2C(x, y, z)dz 0 C(x, y,z)dz 14) 进而源强的积分公式可以根据假设得出： QuC(x, y, z)dydz 15) 16) 17) Q C(x, y,z) exp y2z2 18) 把（ 13）式代入（ 14）积分可以得出： 1 a 2 y 2 1 b b 2 z2 将（13）式和（ 16）式代入（ 15）式可以得出： Q A(x) 最后再将把（ 16）、（17）式代入（ 13）式可以得出： 上式为无界空间连续点源扩散的高斯模型， 然而在实际中， 由于地面的影响， 烟羽是有界的。 根据假设可

20、以把地面看做一镜面， 对泄漏的气体起反射作用， 同 样我们可以利用“像源法”进行处理，原理和示意图在模型一的修正中提到，因 此我们得出： 实源的贡献为： C1(x,y,z) 2 u y z exp( 22 1 y21 (z H )2 2 y2 )exp( 2 2 ) 2 z2 19) 像源的贡献为： Q C2(x, y, z) 2 u y zexp( 1 y2 2 y2 )exp( 2 z2 1 (z H)2) 20) 21) 则该处的实际浓度为： C (x, y, z) C1(x, y, z) C 2( x, y, z) 综合上面的公式得到连续点源高斯烟羽扩散模型： C(x, y, z,t,

21、H ) 2 y exp(2 )exp 2y 2 (z H )2 exp (z H )2 2 z2 22) 5.2.2 模型二的修正 连续点源高斯扩散模型虽然能分析风力对浓度的影响情况， 但为了能更准确 的探究风力对核电站周边放射性物质浓度的影响情况， 我们将考虑泄漏源有效高 度、放射性物质自身的沉降作用和雨水吸附作用对放射性物质浓度的影响， 进而 对连续点源高斯扩散模型进行修正。 5.2.2.1 泄漏源有效高度对模型的修正 如图 2 所示，泄漏源的有效高度 h 是由两部分组成，一是核泄漏口距有效地 面的高度 H ；二是在实际核扩散中核泄漏气团从泄漏口排出时，由于受到热力 抬升和本身动力抬升，进

22、而产生的一个附加高度 H 。因而 h H H 。 H 图四： 烟云抬升示意图 附加高度 H ，主要由核泄漏处泄漏气体的气流具有一初始动量(使他们继 续垂直上升) 和气流温度高于环境温度产生的静浮力决定， 这两种动力引起的烟 云浮力运动称烟云抬升，附加高度 H 即烟云抬升高度，烟云抬升有利于降低地 面的污染物浓度。而且 H 还受到风速、地形地貌等多种因素的影响。 A：当大气稳定度为中性时，计算烟气抬升高度利用 Holland 公式 H vsD (1.5 2.7 Ts Ta D) 1 (1.5vsD 9.6 10 3QH )(23) uTsu 式中： u ：泄漏源出口处的平均风速， m/s； vs

23、：放射性气体出口流速，已知为 mm/s； D ：泄漏源出口的有效内径； QH ：泄漏源的热排放率， Ts：泄漏源出口处温度， Ta ：环境大气平均温度， B：当大气条件为不稳定时，利用 当 QH 21000kw时 x 10H x 10H 当 QH 21000kw时 x 3x* x 3x* kw ； K； K ，取当地近五年的平均值； Briggs 公式计算烟气抬升高度 H 0. 3 6Q2H1 / 3 x 2/ 3 u( ) H 1. 5Q5H1 / 3 H 2 / 3 (u ) H 0. 3 6Q2H1 / 3 x1 / 3 u( ) H 0. 3 3Q2H 3 / 5 H 2/ x* 0.

24、33QH 3/5 H 3/5 (u) 6/5 综上所述，泄漏源的有效高度为： 24) gd2 25) 5.2.2.2 考虑放射性物质自身的沉降作用对模型的修正 放射性物质的沉降速度取决于空气阻力和自身重力， 利用斯托克斯公式表示 沉降速度： Vs 18 ：放射性物质粒子密度， kg/ m3； g ：重力加速度， 9. 806 5m/s ； d ：放射性物质粒子直径， m ； ：空气的动力粘性系数，可取 1.8 10 5 kg /(m s) ； Vs ：沉降速度， m/ s ，含碘放射性核素的沉降速度为 Vs= 1. 1 cm/ s； 在扩散过程中重力沉降的位移叠加在羽流中心线上， 使中心线向下

25、倾斜， 放射性 物质粒子则相当于在下倾的中心线上扩散， 放射性物质的扩散与沉降的叠加可认 为是放射源以 Vs的速度向下移动。在 x处向下移动的垂直距离为 Vst Vsx ，即泄 s 漏源的有效高度 h下降了 Vst Vsx，泄漏源的有效高度成为 h Vst h Vsx ，考 虑到地面的全反射作用，反射项的有效高度也变成了 h Vst h Vsx 。则修正后 的连续点源高斯扩散模型为： C(x, y,z,h) exp( 2 y2 )exp ( z h ( z h exp 5.2.2.3 考虑雨水的吸附作用对模型的修正 降雨对放射性物质的浓度有一定影响， 即雨水对放射性物质有一定的吸附作 用。以吸

26、附系数 来表示雨水对放射性物质吸附作用的大小， 与降雨强度的关 系为：aI b，式中 I 为降雨强度， a,b为经验系数。如果放射性物质含碘，则 a 8 10 5,b 0.6；反之， a 1.2 10 5 ,b 0.5。 雨水的吸附作用导致的放射性物质浓度的减小，可对源强进行修正： x Q(x) Q exp()(27) 则进一步修正的连续点源高斯扩散模型为： ( z h Vsx)2 2 z2 Vs x 2 Q( x)y2(z h s )2 C(x, y,z,h) exp(2 )exp2 exp 2 y22 y2 28) 综上所述：修正的连续点源高斯扩散模型为： exp (z h Vsx )2

27、(z h2 z2 ) Vsx 2 Q(x)y2( z h s )2 C(x, y, z, h)exp( 2 )exp 2 2 u y z 2 y22 y2 Q(x) Q exp( x),aI b Vs gd2 18 29) 5.2.3 模型二的求解 模型所需参数的选取对模型的求解至关重要， 通常情况下气象参数的选取是 利用该地区多年气象资料，采取工业安全与环保统计的方法进行有关参数的确 定，而其他扩散参数则以实际测定为准。 A：大气稳定度的计算 根据我国标准 (GB/ T 1320191) )制订地方大气污染物排放标准的技术方法 的规定，大气稳定度分为 6级，分别为 A极不稳定、 B不稳定、

28、C弱不稳定、 D中性、 E弱稳定、 F稳定。该方法的技术路线是：根据核泄漏源所在地的 经度和纬度以及泄漏的日期和时间计算当时的太阳高度角 h0 ，利用天气条件确定 辐射等级，然后利用辐射等级和风速确定大气稳定度， 最后查扩散参数幂函数表， 确定扩散参数。 首先，然后，由太阳高度角 h0 和云量查出太阳辐射等级；最后，再根据地面 风速确定当时的大气稳定度。 总云量/ 低云量 夜间 太阳高度角 h0 h0 15 15 h0 35 35 h0 65 h0 65 4/ 4 -2 -1 +1 +2 +3 5 7/ 4 -1 0 +1 +2 +3 8/ 4 -1 0 0 +1 +1 7/57 0 0 0

29、0 +1 8/ 8 0 0 0 0 0 表一：太阳辐射等级 地面风 速 (m/s) 太阳辐射等级 +3 +2 +1 0 -1 -2 1.9 A A-B B D E E 2 2.9 A-B B C D E F 3 4.9 B B-C C D D E 5 5.9 C C-D D D D D 6 C D D D D D 注：地面风速系指距地面 10m 高度处 10min 平均风速 表二：大气稳定度的等级 B：扩散参数的确定 扩散参数 y, z 的确定，采用 Briggs 给出一套扩散参数幂函数表，如表三和 表四： 大气稳定度 y z A 0.22x(1 0.0001x) 2 0.2x B 0.16x

30、(1 0.0001x) 2 0.12x C 0.11x(1 0.0001x) 12 0.08x(1 0.0002x)12 D 0.08x(1 0.0001x) 12 0.06x(1 0.0015x ) 1 2 E 0.06x(1 0.0001x) 2 0.03x(1 0.0003x) F 0.04x(1 0.0001x) 2 0.016x(1 0.0003x) 表三： Briggs 扩散参数(平原地区和城市远郊区) 大气稳定度 x y A-B 1 0.32x(1 0.0004x) 2 1 0.24x(1 0.0001x) 2 C 0.22x(1 0.0004x) 2 0.20 x D 0.16

31、x(1 0.0004x) 2 0.14x(1 0.0003x) 2 E-F 0.11x(1 0.004x) 12 0.08x(1 0.0015x) 12 表四： Briggs 扩散参数(工业区和城市中心区) 5.3 模型三的建立与求解 当风速为 km/ s时，建立上风和下风 L公里处放射性物质浓度的预测模型。 此问题可以利用模型二进行求解，将模型二中的 u 替换成 k s,其中 k 为风速的 向量表示， s为泄漏点放射性气体自身扩散速度的向量表示。 根据假设，风速 k 沿 x轴正方向，恒为正。对上风向 L 公里处放射性气体浓度计算时， s方向沿 x轴 负方向，为负；对下风向 L 公里处放射性气

32、体浓度计算时， s方向沿 x 轴正方向， 为正。因此此问题转化为求 (L,0, z) 和( L,0, z)两处的放射性气体浓度， 即风速为 km/ s时，上风和下风 L 公里处放射性物质浓度的预测模型为： Q(x) 2 ( k s) y z C( L,0,z) = exp( 2 y 2y )exp 2 (kQ+(sx)exp( 2 (k +s) y z 2 y 2y )exp Vsx 2Vsx 2 ( z h s )2(z h s )2 2 exp 2, 上风向 L处 2 z 2 z (z h Vsx)2(z h Vsx)2 2 exp 2，下风向 L处 2 z 2 z 30) 我们根据模型分

33、别模拟下风向和上风向的放射性物质扩散仿真图，如图 五、图六：（程序分别见附录二、三） 图五：核泄漏下风向扩散图 北京： J 内 环境 图六：核泄漏上风向扩散图 七 参考文献 1 卓金武，魏永生，秦建，李必文 Matlab 在数学建模中的应用 M 北京航空航天大学出版社， 2011 2 程勇，于林，姚安林采用高斯模型分析输气管道泄漏后气体的扩散 蒙古石油化工， 2010 ，14： 49-51 3 张斌才，赵军大气污染扩散的高斯烟羽模型及其 GIS 集成研究 J 监测管理与技术， 2008， 20（5）：17-19 八 附录 附录一： %绘制当没有风影响时放射性物质的扩散图形 function =expand(k,Q,C,t,r,ws,wa,wd) R=sqrt(-4)*k*t*log(C/(1-r)*Q).*(4*pi*k*t)(1.5) for i=0:24:t/3600 if (wd*(i/24)=R R=R-wd; R0=R; t=0:0.1:0.5*p