重金属污染分析扩散模型_XXXX建模.docx

2011数学建模扩散模型承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从D/E/F中选择一项填写）： 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)： 日期： 年 月 日赛区评阅编号（由赛区组委会评阅前进行编号）：2011重庆邮电大学数学建模第二次模拟编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：放射性气体扩散浓度预估模型【摘要】本文是以日本地震引起的福岛核电站的核泄漏为背景，并以给出的数据为基础，研究某一假设核电站的核泄漏问题。我们通过收集相关的资料，并结合题目给出的数据，建立了高斯模型、连续点源高斯扩散模型解决了题目提出的四个问题。针对问题一：考虑到泄漏源是连续、均匀和稳定的，我们运用散度、梯度、流量等数学概念，通过“泄漏放射性物质质量守恒”、“气体泄漏连续性定理”、 Guass公式及积分中值定理得到了无界区域的抛物线型偏微分方程，然后再通过电源函数解出空间任意一点的放射性物质浓度的表达式，把此表达式定为模型一的前身。鉴于放射性物质的扩散受到诸多因素的影响，如：泄漏源的实际高度、地面反射等。我们以泄漏口为坐标原点建立三维坐标系，通过“像源法”处理地面反射对放射性物质浓度的影响，并由此对模型一的前身进行修正完善，得到模型一：高斯模型，即放射性物质浓度的预测模型。最后我们模拟了放射性物质无风扩散仿真图。针对问题二：当风速为k m/s时，我们根据放射性核素云团在大气中迁移和扩散的数值计算的基本方法和步骤，并以泄漏点源在地面的投影点为坐标原点，以风向方向为轴，铅直方向为轴，与轴水平面垂直方向为轴建立三维坐标系，地面的反射作用同样利用“像源法”进行处理，得到连续点源高斯扩散模型。考虑到地面反射、烟云抬升、放射性物质自身的沉降及雨水的吸附等对浓度的影响，我们对连续点源高斯扩散模型进行了修正，建立了修正的连续点源高斯扩散模型。最后利用大气稳定度确定了扩散参数，进而求解了模型。针对问题三：经分析，问题三的提出是以问题二为基础的，模型三的建立只需要将模型二加以调整即可。我们以风速方向为轴正方向，将风速与放射性物质的扩散速度进行矢量运算，此问题则转化为求和两点处的放射性物质浓度，由此建立模型三，即上风和下风公里处放射性物质浓度浓度的预测模型。针对问题四：首先，我们通过网络收集了相关数据，然后，我们结合模型二、模型三对数据进行整理代入，算出了日本福岛核电站泄漏的放射性物质扩散到中国东海岸和美国西海岸的浓度分别为、。关键词：高斯模型 连续点源高斯扩散模型 核泄漏 一 问题的重述1.1 问题背景目前，核电站的发展能带来巨大的经济效益和社会效益，但核电站一旦发生核泄漏，将会给人们的生命健康和周边环境带来巨大的危害性影响。2011年3月日本的福岛核电站的放射性气体的核泄漏事件更让我们关注放射性气体泄漏时的浓度问题。因此，正确的测出大气中放射性物质的浓度在环境监测和安全评估中具有重要意义1.2 问题提出有一座核电站遇自然灾害发生泄漏，浓度为的放射性气体以速度m kg/s匀速排出，在无风的情况下，以速度s m/s匀速在大气中向四周扩散.1） 在无风的情况下，建立一个描述核电站周边不同距离地区、不同时段放射性物质浓度的预测模型。2） 当风速为k m/s时，给出核电站周边放射性物质浓度的变化情况。3） 当风速为k m/s时，分别给出上风和下风L公里处，放射性物质浓度的预测模型。4） 将建立的模型应用于福岛核电站的泄漏，计算出福岛核电站的泄漏对我国东海岸，及美国西海岸的影响。计算所用数据可以在网上搜索或根据具体情况自己模拟。二 符号说明空间任意一点的放射性物质的扩散系数空间任意一点的放射性物质浓度放射性气体的扩散速度, m/ s泄漏源泄漏的放射性物质总量空间域空间域其体积一规则的球面面积在内通过的流量内放射性物质的增量任意扩散时刻泄漏源距地面的实际高度烟云抬升高度泄漏源有效高度泄漏源高度处的平均风速，m/ s源强，用浓度标准差表示的轴上的扩散参数放射性气体出口流速,m/s泄漏源出口的有效内径泄漏源的热排放率，泄漏源出口处温度，环境大气平均温度，沉降速度，空气的动力粘性系数雨水吸附系数太阳高度角三 模型假设1、扩散过程中浓度在、轴上的变化分布是高斯分布。2、放射性物质的扩散看作是空间某一连续点源向四周等强度地瞬时释放放射性物质，放射性物质在无穷空间扩散过程中不发生性质变化，且不计地形影响。3、放射性物质扩散服从扩散定律，即单位时间通过单位法向面积的流量与它的浓度梯度成正比。4、放射性物质在穿过降雨区域时，其强度由于雨水的吸收而减少，减少比率为常数。5、假设地面对放射性气体起全反射作用，地面和海面对放射性气体没有吸附，将海面视为平原地区6、假设风向为水平风向，且风向风速不随时间变化。7、扩散过程中不考虑泄漏点内部温度的变化对气体扩散的影响。 四 问题的分析4.1 问题（1）的分析一座核电站遇自然灾害发生核泄漏，浓度为放射性气体以速度m kg/s匀速排除，这近似于放射性物质源是连续均匀稳定的。在无风情况下，放射性气体以速度s m/s匀速在大气中向四周扩散，放射性气体的扩散服从扩散定律，即单位时间通过单位法向面积的流量与它的浓度梯度成正比。在这些条件下，我们明确了要研究的问题是点源连续泄漏的扩散问题，题给条件中明确要求不考虑风力的影响，但为了使建立的模型更加贴近实际，需考虑地面反射、核泄漏源的实际高度、降雨等因素对浓度分布的影响。由“扩散定律”“放射性物质质量守恒定律”“气体泄漏连续性定理”可得出无界区域的抛物线型偏微分方程。再通过假设条件建立未考虑地面反射、核泄漏源的实际高度、降雨等浓度影响因素的初步模型，然后从这些影响因素对模型进行完善，最终得出核电站周边不同距离地区、不同时段放射性物质浓度的预测模型。4.2 问题（2）的分析本问是探究风速为时，核电站周边放射性物质浓度的变化情况。当环境中空气流动产生风力时，在均匀湍流场中，扩散参数与下风向距离的关系是明确的，核泄漏泄漏时间较长时，可认为扩散是稳定的。在下风向的湍流扩散相对于风力引起的移流相可忽略不计，在流动方向建立x轴，不考虑横向速度和垂直速度。根据假设，空间中放射性物质的浓度服从高斯分布，可利用连续点源放射性物质的高斯扩散模型。放射性物质在大气中扩散受诸多因素影响，考虑泄漏源有效高度、放射性物质自身重力产生的重力沉降、雨水的吸附等因素对放射性物质浓度的影响是必要的，通过这些影响因素对高斯模型进行修正，然后利用修正后的高斯模型探究核电站周边放射性物质浓度的变化情况。4.3 问题（3）的分析本问是要求当风速为时，建立上风和下风公里处的放射性物质浓度的预测模型。经分析，此问是问题（2）的延伸，我们只需建立合适的坐标系，将此问题转化为求具体两处的放射性物质的浓度，便能得出上风和下风公里处的放射性物质浓度预测模型。4.4 问题（4）的分析本问要求将之前建立的模型应用于日本福岛核电站的泄泄漏，计算出福岛核电站的泄漏中国东海岸和美国西海岸的影响，此问实际上是将模型二、模型三具体化。通过收集求解模型的相关数据，利用模型和数据模拟日本福岛核泄露对中国东海岸和美国西海岸的影响。五 模型的建立与求解5.1 模型一的建立与求解5.1.1 模型一的初步建立以核泄漏点正下方的地面为坐标原点，平均风向为X轴、指向下风方向，铅直方向为Z轴，水平垂直于风向轴（X轴）为Y向，建立空间坐标系，则核电站泄漏点距有效地面的高度为，则泄漏点位置坐标为。 图一：空间坐标系示意图将气体从泄漏源泄漏时刻记作t=0，时刻t无穷空间中任意一点坐标为（x,y,z）的浓度记为C（x,y,z,t）,根据假设2，单位时间通过单位法向面积的流量与浓度梯度成正比，则： （1）是扩散系数，表示梯度，负号表示由浓度高向浓度低的地方扩散。考察空间域，其体积为，包围的曲面为，为一规则的球面，外法线向量为。则在内通过的流量为： （2）内放射性物质的增量为： （3）从泄漏源泄漏的放射性物质的总量为： （4）根据“质量守恒定律”和“气体泄漏连续性原理”，单位时间内通过所选曲面的向外扩散的放射性物质的量与曲面内放射性物质增量之和，等于泄漏源在单位时间内向外泄漏的放射性物质。则： （5）即：（6）又根据曲面积分的Gauss公式： (其中是散度记号) （7）由以上两式得：即为： （8）由以上公式并利用积分中值定理得： （9） 这是无界区域的抛物线型偏微分方程，根据假设1，初始条件为作用在坐标原点的电源函数，记作 （10） 方程（9）满足方程（10）的解为： 此模型只是在不考虑风速的情况下建立的，但为了使模型具有更加的实用性，下面我们将考虑泄漏源的实际高度、地面反射、降雨等因素对浓度的影响，完善模型。 5.1.2 模型一的最终建立1. 地面反射对模型的完善泄漏源有一定的高度，且泄漏点源是连续点源，则泄漏点源可视为高架连续点源，考虑到地面对扩散来的放射性气体有反射作用，根据假设4，地面对到达地面的扩散气体完全反射。这儿可认为地面就像镜子一样，对放射性气体起全放射作用，可用“像源法”处理，如图3，建立三个坐标系，一是以泄漏源（实源）为坐标原点；二是以泄漏源在地面的投影点为原点，点是空间的任意一点，坐标为；三是以泄漏源关于地面的像对称源（像源）为原点。把点放射性气体浓度看成两部分（实源与像源）作用之和。图二： 高架连续点源扩散示意图从以上分析知，点放射性气体的浓度为实源和像源的放射性气体扩散至此点浓度的叠加。则实际泄漏源（实源）对点的浓度贡献部分可用来表示；因为地面对扩散物质完全反射，则像对称源（像源）对点的浓度贡献部分可用来表示。于是对（11）式所修正完善的模型为： （12）我们自己模拟一组数据，利用matlab进行仿真模拟，可实现该模型的模拟图像。假设扩散系数=0.00001，放射性物质的初始浓度=100，扩散时间t=1000000，放射源总量=1000000。（程序见附录一）图三：核泄漏无风扩散5.2 模型二的建立与求解5.2.1 模型二的建立放射性核素云团在大气中迁移和扩散的数值计算基本上可分为二步。第一步根据大气动力学理论进行所关心区域中风场的计算，其理论基础是大气运动方程、连续性方程、状态方程、热力学方程和水汽方程构成的基本方程组。在大气科学研究领域中，已有多个实用的大气环流模式。第二步进行已知风场中放射性核素云团迁移和扩散的计算，可采用类似于处理大气污染的方法，假设放射性核素云团不影响大气流体速度和温度，求解放射性核素的连续性方程。当风速为时，利用连续点源高斯扩散模型分析核电站周边放射性物质浓度的变化情况。此泄漏点源是有边界点源，泄漏点源的实际高度为。以泄漏点源在地面的投影点为坐标原点，以风向方向为轴，铅直方向为轴，与轴水平面垂直方向为轴建立三维坐标系，由于扩散过程中浓度在、轴上的变化分布符合高斯分布，所以下风向的任意一点的浓度函数为： （13）根据概率统计我们可以得出方差的表达式为： （14）进而源强的积分公式可以根据假设得出： （15）把（13）式代入（14）积分可以得出： （16）将(13)式和（16）式代入（15）式可以得出： （17） 最后再将把（16）、（17）式代入（13）式可以得出： （18）上式为无界空间连续点源扩散的高斯模型，然而在实际中，由于地面的影响，烟羽是有界的。根据假设可以把地面看做一镜面，对泄漏的气体起反射作用，同样我们可以利用“像源法”进行处理，原理和示意图在模型一的修正中提到，因此我们得出：实源的贡献为： （19）像源的贡献为： （20）则该处的实际浓度为： （21）综合上面的公式得到连续点源高斯烟羽扩散模型： （22）5.2.2 模型二的修正连续点源高斯扩散模型虽然能分析风力对浓度的影响情况，但为了能更准确的探究风力对核电站周边放射性物质浓度的影响情况，我们将考虑泄漏源有效高度、放射性物质自身的沉降作用和雨水吸附作用对放射性物质浓度的影响，进而对连续点源高斯扩散模型进行修正。5.2.2.1 泄漏源有效高度对模型的修正如图2所示，泄漏源的有效高度是由两部分组成，一是核泄漏口距有效地面的高度；二是在实际核扩散中核泄漏气团从泄漏口排出时，由于受到热力抬升和本身动力抬升，进而产生的一个附加高度。因而。图四： 烟云抬升示意图附加高度，主要由核泄漏处泄漏气体的气流具有一初始动量（使他们继续垂直上升）和气流温度高于环境温度产生的静浮力决定，这两种动力引起的烟云浮力运动称烟云抬升，附加高度即烟云抬升高度，烟云抬升有利于降低地面的污染物浓度。而且还受到风速、地形地貌等多种因素的影响。A：当大气稳定度为中性时，计算烟气抬升高度利用公式 （23）式中： ：泄漏源出口处的平均风速，m/s；：放射性气体出口流速，已知为m/s；：泄漏源出口的有效内径；：泄漏源的热排放率，；：泄漏源出口处温度，；：环境大气平均温度，取当地近五年的平均值；B：当大气条件为不稳定时，利用公式计算烟气抬升高度当时 当时 综上所述，泄漏源的有效高度为： （24） 5.2.2.2 考虑放射性物质自身的沉降作用对模型的修正放射性物质的沉降速度取决于空气阻力和自身重力，利用斯托克斯公式表示沉降速度： （25）：放射性物质粒子密度，；：重力加速度，9. 806 5；：放射性物质粒子直径，；：空气的动力粘性系数，可取；：沉降速度，含碘放射性核素的沉降速度为= 1. 1 cm/ s；在扩散过程中重力沉降的位移叠加在羽流中心线上，使中心线向下倾斜，放射性物质粒子则相当于在下倾的中心线上扩散，放射性物质的扩散与沉降的叠加可认为是放射源以的速度向下移动。在处向下移动的垂直距离为，即泄漏源的有效高度下降了，泄漏源的有效高度成为，考虑到地面的全反射作用，反射项的有效高度也变成了。则修正后的连续点源高斯扩散模型为：5.2.2.3 考虑雨水的吸附作用对模型的修正降雨对放射性物质的浓度有一定影响，即雨水对放射性物质有一定的吸附作用。以吸附系数来表示雨水对放射性物质吸附作用的大小，与降雨强度的关系为：，式中为降雨强度，为经验系数。如果放射性物质含碘，则；反之，。雨水的吸附作用导致的放射性物质浓度的减小，可对源强进行修正： （27）则进一步修正的连续点源高斯扩散模型为：（28）综上所述：修正的连续点源高斯扩散模型为：（29）5.2.3 模型二的求解模型所需参数的选取对模型的求解至关重要，通常情况下气象参数的选取是利用该地区多年气象资料，采取工业安全与环保统计的方法进行有关参数的确定，而其他扩散参数则以实际测定为准。A：大气稳定度的计算根据我国标准(GB/ T 1320191) )制订地方大气污染物排放标准的技术方法的规定，大气稳定度分为6级，分别为A极不稳定、B不稳定、C弱不稳定、D中性、E弱稳定、F稳定。该方法的技术路线是：根据核泄漏源所在地的经度和纬度以及泄漏的日期和时间计算当时的太阳高度角，利用天气条件确定辐射等级，然后利用辐射等级和风速确定大气稳定度，最后查扩散参数幂函数表，确定扩散参数。 首先，然后，由太阳高度角和云量查出太阳辐射等级；最后，再根据地面风速确定当时的大气稳定度。总云量/低云量夜间太阳高度角-2-1+1+2+3-10+1+2+3-100+1+10000+100000表一：太阳辐射等级地面风速太阳辐射等级+3+2+10-1-21.9AA-BBDEEA-BBCDEFBB-CCDDECC-DDDDD6CDDDDD注：地面风速系指距地面10m高度处10min平均风速表二：大气稳定度的等级B：扩散参数的确定扩散参数的确定，采用给出一套扩散参数幂函数表，如表三和表四：大气稳定度ABCDEF表三： 扩散参数（平原地区和城市远郊区）大气稳定度A-BCDE-F表四： 扩散参数（工业区和城市中心区）5.3 模型三的建立与求解当风速为时，建立上风和下风公里处放射性物质浓度的预测模型。此问题可以利用模型二进行求解，将模型二中的替换成,其中为风速的向量表示，为泄漏点放射性气体自身扩散速度的向量表示。根据假设，风速沿轴正方向，恒为正。对上风向公里处放射性气体浓度计算时，方向沿轴负方向，为负；对下风向公里处放射性气体浓度计算时，方向沿轴正方向，为正。因此此问题转化为求和两处的放射性气体浓度，即风速为时，上风和下风公里处放射性物质浓度的预测模型为： （30）我们根据模型分别模拟下风向和上风向的放射性物质扩散仿真图，如图五、图六：（程序分别见附录二、三）图五：核泄漏下风向扩散图图六：核泄漏上风向扩散图5.4 问题（4）的求解此问题可结合模型二和模型三求解，我们通过网络收集了相关数据，将这些数据代入模型进行模拟。数据显示，日本在核泄漏时期的风向主要为偏西风，风速为2-3级,辐射源强度。建立如图七（东海岸扩散图）中所示的坐标轴，则中国东海岸的坐标大概可取为（-1500km,1900km）。建立如图八（西海岸扩散图）中所示的坐标轴，则美国西海岸坐标大概取为（4300,7500）。根据模型二的求解方法可确定，风速为4m/s，分别求解福岛核电站核泄漏泄漏的放射性物质扩散到中国东海岸和美国西海岸的浓度：中国东海岸：=美国西海岸：=图七：东海岸扩散图图八：西海岸扩散图六 模型的评价及改进6.1模型的评价6.1.1 模型优点第一：模型一，我们考虑了泄漏源高度、地面反射等因素对放射性物质浓度的影响，利用这些影响因素对浓度预测模型进行了修正完善，建立的浓度预测模型具有实用性、高精度性、合理性。第二：模型二，我们建立了连续点源高斯扩散模型，并考虑了烟云抬升高度，求得泄漏源的有效高度，通过泄漏源的有效高度、放射性物质自身的沉降作用、雨水的吸附作用对模型进行修正，得到了更优化的连续点源高斯扩散模型。此模型能较好的求出核电站周边放射性物质浓度，并与模型一求出的浓度进行比较，从而给出核电站周边放射性物质浓度的变化情况。第三：模型三是在模型二的基础上加以修改转化的，因此模型三继承了模型二的优点，使放射性物质浓度的预测更加精确。第四：我们建立的以上模型在核事故的应急救援过程中，对救援人员划定警戒区和确定周围居民的疏散范围具有重要意义，并可为制定救援方案和应急决策提供科学依据。对放射性云团在空中迁移和扩散提供浓度的定量描述，也对可能发生的核事件的放射性核素浓度监测及监测时间范围提供相关信息。6.1.2 模型缺点第一：模型只考虑了平原和城市两种情况，不利于扩散参数的确定。但地形复杂多变，不能找到充分的数据对其进行分析。第二：模型缺乏检验。没有相关数据和标准用于模型的检验。6.2 模型的改进日本是多山国家，考虑山林对放射性物质扩散的影响能更好的预测和估算放射性物质的浓度；放射性物质的扩散，特别是向其他国家地区扩散，放射性物质的扩散大多要跨洋越海，考虑海水对放射性物质的吸收显得至关重要，将此因素考虑进去，能使模型更加贴近实际。七 参考文献1 卓金武，魏永生，秦建，李必文Matlab在数学建模中的应用 M 北京：北京航空航天大学出版社，20112 程勇，于林，姚安林采用高斯模型分析输气管道泄漏后气体的扩散 J 内蒙古石油化工，2010 ，14：49-513 张斌才，赵军大气污染扩散的高斯烟羽模型及其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/3600if (wd*(i/24)<=R R=R-wd; R0=R;t=0:0.1:0.5*pi;p=0:0.1:2*pi;wsx=ws*cos(wa)*i;wsy=ws*sin(wa)*i;theta,phi=meshgrid(t,p);x=R*sin(theta).*cos(phi)+wsx;y=R*sin(theta).*sin(phi)+wsy;z=R*cos(theta);hold onsurfc(x,y,z) endendR0=R0%绘制其它相关信息view(20,70) %观察角度title('核泄漏无风影响扩散');%图形标题xlabel('x'); %x轴标记ylabel('y'); %y轴标记zlabel('z'); %z轴标记axis(-10,15,-10,20,0,12) %xyz轴显示范围余量x0=ws*cos(wa)*1000000/3600y0=ws*sin(wa)*1000000/3600R1=sqrt(x0)2+(y0)2)End附录二：Q=2; %输入源强k=2; %输入风速s=0.5;u=k+s;d=1; %步长Z=0.4; x=10:d:300; %下风向距离y=-70:d:70; %横风向距离x,y=meshgrid(x,y);by0=0.08*x.*(1+0.0001*x).(-1/2);bz0=0.06*x.*(1+0.0015*x).(-1/2);tempy1=-y.*y./by0./by0./2;tempy2=2.718282.(tempy1);c=Q/pi/u*(by0.*bz0).(-1).*tempy2Cs=2000; %输入求解的条数contour(x,y,c,Cs);shading interp;colorbar;grid;title('核泄漏模拟图一')附录三：Q=2; %输入源强k=2; %输入风速s=0.5;u=k-s;d=1; %步长Z=0.4; x=-300:d:-10; %下风向距离y=-70:d:70; %横风向距离x,y=meshgrid(x,y);by0=0.08*x.*(1+0.0001*x).(-1/2);bz0=0.06*x.*(1+0.0015*x).(-1/2);tempy1=-y.*y./by0./by0./2;tempy2=2.718282.(tempy1);c=Q/pi/u*(by0.*bz0).(-1).*tempy2Cs=2000; %输入求解的条数contour(x,y,c,Cs);shading interp;colorbar;grid;title('核泄漏模拟图二')