数学建模优秀论文-放射性气体扩散预测.doc

摘 要本文是以日本福岛核电站遭遇自然灾害发生核泄漏的背景而提出的。对于问题一，考虑到放射性物质的泄漏是连续不断的。本文根据“泄漏放射性物质质量守恒定律”和“气体泄漏连续性原理”建立了微积分方程，应用了高等数学中散度、梯度、流量等数学概念，通过Guass公式、四维二阶偏微分方程，因而得到了核电站周边不同距离地区、不同时段放射性物质浓度的预测模型。同时为使模型适用范围更广，本文引入了地面反射系数，考虑了由于放射性物质从泄漏口喷出时具有初动量而使其泄漏源有效高度被抬高等因素，进而得到了在无风环境中适用范围更广的“高斯修正模型”。对于问题二，要探究风速对放射性物质浓度分布的影响。本文运用概率学1知识，通过图解和数学推导得出“连续点源放射性物质高斯扩散模型”。本文依次考虑了“重力沉积”、“雨水沉积”、“核衰变”等因素对浓度分布的影响。并通过构建“耗减因子”、“衰变因子”等方法将耗减和衰变的放射性物质“投影”到泄漏源浓度中，得到了经多次合理修正后的“优化高斯模型”，并据此分析了泄漏源周边地区放射性物质的浓度变化。针对于问题三，本文在问题二的基础上，结合考虑风速和放射性物质扩散速度在空间中的矢量运算。得出在对上风口分析时，要分类讨论风速和自然扩散速度之间的大小关系，当风速小于自然扩散速度时，放射性物质是无法到达上风口的。对于问题四，本文参阅整理大量气象、地理、新闻资料，选择我国东海岸典型地域-山东半岛作为研究对象，综合考虑对应海域平均风速及风向、地理距离、海水对放射性物质扩散的部分反射系数等因素，并通过C+编程模拟计算，预测出放射性核物质将经过6.5天到达我国东海岸，且131I浓度预测值为：0.100与实际情况比较吻合。关键词：放射性气体 扩散 浓度变化 高斯修正模型 预测 1 问题的提出由于重大的突发性核泄漏紧急灾害事件具有爆发性、空间分布不连续性、对周边地形和气象条件的敏感性的特点，研究核事故所释放的物质的时空分布需要高度精确的技术，但是在对于更好地保护环境有着极其重要的意义。在有一座核电站遇自然灾害发生泄漏，浓度为的放射性气体以匀速排出，速度为，在无风的情况下，匀速在大气中向四周扩散，速度为。 问题一，若能建立一个描述核电站周边不同距离地区、不同时段放射性物质浓度的预测模型，这对于研究核污染模式具有重要的意义。 问题二，当风速为时，给出核电站周边放射性物质浓度的变化情况，这对于研究核电站附近浓度的在实际环境下有着重要的作用。 问题三，当风速为时，计算出上风和下风公里处的放射性物质浓度的预测模型就显得更加急迫的了。 问题四，将建立的模型应用于福岛核电站的泄漏，计算出福岛核电站的泄漏对我国东海岸的影响，这一实际的有意义预测，可以明确我们实际的污染情况，为我们的核应急决策提供技术支持。2 问题的分析对于问题一，在无风的情况下，放射性气体以m/s的速度，匀速在大气中向四周扩散。在此条件下，探求一个模型来对核电站周边不同距离地区、不同时段放射性物质浓度进行预测。我们要明确此问题研究的核扩散是点源连续泄露的扩散问题。虽然只是要求考虑在无风情况下放射性物质浓度分布，但为了使模型更贴切实际，需考虑地面反射、泄漏源有效高度等因素对浓度分布的影响。根据“泄露放射性物质质量守恒定律”和“气体泄漏连续性原理”进行分析，发现要得出核电站周边不同距离地区、不同时段放射性物质浓度的预测模型，最后对于该方程进行分析求解。对于问题二，为了探究风速对发生核泄漏的核电站周边放射性物质浓度分布的影响，运用概率学知识，通过图解和数学推导得出“连续点源放射性物质高斯扩散模型”。应在“连续点源放射性物质高斯扩散模型”的基础上经多次合理修正后得到更好的“优化高斯模型”。对于问题三，该问题要求建立泄漏源上风口和下风口处放射性物质浓度的预测模型，在参考第二问的基础上，主要考虑风速和放射性物质扩散速度在空间中的矢量运算。在对上风口分析时，要分类讨论风速和自然扩散速度之间的大小关系，当风速小于自然扩散速度时，放射性物质是无法到达上风口的。对于问题四，应参考大量气象、地理、新闻资料，选择我国东海岸典型地域-山东半岛，作为研究对象，综合考虑对应海域平均风速及风向、地理距离、海水对放射性物质扩散的部分反射系数等因素，预测出放射性核物质与实际情况比较。3 模型的假设考虑到放射性气体扩散的复杂性，为简单起见，在讨论扩散模型时都作了如下假设:（1）瞬时泄漏假定瞬时完成，连续泄漏假定泄漏速率恒定；（2）气云在平整、无障碍物的地面上空扩散；（3）气云中不发生化学反应，地面对气云无吸收；（4）为水平风向，风速和风向不随时间变化。4 符号说明及名词解释4.1符号说明符号说明放射性气体的传播速度风速，单位泄漏点距有效地面的高度任意扩散时刻空间任意一点的放射性物质浓度空间任意一点的放射性物质的扩散系数空间域空间域其体积一规则的球面面积在内通过的流量内放射性物质的增量符号说明从泄漏源泄漏的放射性物质的总量附加高度核泄漏出口处的温度环境温度设地面反射系数源强，单位为分别为用浓度标准差表示的轴上的扩散参数沉降速度，单位为地面干沉积率冲洗系数放射性核素的半衰期4.2名词解释烟羽又称烟云（smoke cloud)、烟流（smoke plume)：从烟囱中连续排放到大气中的烟气流。由于烟羽各部分的运动速度不同，因而其外形也千变万化。不同的烟羽形状表示污染物浓度的空间分布不同。它与大气湍流、大气稳定度、地形地物、排放参数等有密切的关系。动力抬升：暖气流受锋面、辐合气流的作用被迫上抬，或者在运行中受地形阻挡产生上升运动，这种空气在运动中由外力（不包括重力和浮力）使一部分空气被抬上升。湍流扩散：是指湍流运动导致大气或水体中的污染物质或其他物质与周围洁 流体的混合。5 模型的建立和求解5.1问题一模型的建立与求解5.1.1模型一的建立以核泄漏点正下方的地面为坐标原点，平均风向为X轴、指向下风方向，铅直方向为Z轴，水平垂直于风向轴（X轴）为Y向，建立空间坐标系，则核电站泄漏点距有效地面的高度为，则泄漏点位置坐标为。图1 空间坐标系示意图并记时刻时，空间任意一点的放射性物质浓度为。根据假设设单位时间通过单位法向面积的流量2与浓度梯度成正比，有： （1）是扩散系数，表示梯度，负号表示由浓度高向浓度低的地方扩散。先考察空间域，其体积为，包围的曲面为，为一规则的球面，外法线向量为。则在内通过的流量为： （2）内放射性物质的增量为： （3） 从泄漏源泄漏的放射性物质的总量为： （4） 根据质量守恒定律和连续性原理，单位时间内通过所选曲面的向外扩散的放射性物质与曲面内放射性物质增量之和，等于泄漏源在单位时间内向外泄漏的放射性物质。有： （5）即，（6）又根据曲面积分的Gauss公式3： (其中是散度记号) （7）由以上两式得：即为： （8）根据参考文献 （9）解得： （10） 上述模型仅是一个最理想化的预测浓度模型5，因为它将环境视为无边界空间，且没有考虑放射性物质从泄漏口释放出时的初动量。为使所建立模型更加贴合实际，本文先从“有效泄漏源”和“地面反射”这两个方面对所建立模型进行修正。5.1.2考虑热力抬升作用对模型的修正如图2所示，为核泄漏点源的有效源高。它是由两部分构成：一是核泄漏口的有效高度；二是在实际核扩散中核泄漏气团从泄漏口排出时，由于受到热力抬升和本身动力抬升，进而产生的一个附加高度。因而。图2 热力抬升示意图对于，主要由浮升力和泄漏的初始动量决定，同时还要受到泄漏口温度、大气温度、风速、地形地貌等多种因素的影响。我们直接引用气体污染扩散学中应用较广范的，有关烟气抬升高度的综合分析公式6： （11）式中： 泄漏源的实际高度；泄漏源出口处的风速，已知为m/s；泄漏源出口的有效直径；放射性气体的扩散速度，已知为m/s；浮力通量，；由Briggs抬升公式知，浮力通量，其中核泄漏出口处的温度；环境温度。a：在有风()且释放气体温度与环境温度差35K()时，抬升高度： （12）其中，为大气压强b：小风()时，且温度差35K()，抬升高度为： （13）其中为泄漏源的有效高度处上的环境温度梯度()。c：当温度差35K()时，此时的抬升高度： （14）此时，当风速时，取= 结合以上分析，我们将以上结论应用到本题中（带入风速k及放射性气体的扩散速度s）得： 综上所述，泄漏源的有效高度为： （15）5.1.3考虑到地面反射对模型的修正考虑到地面会对扩散来放射性气团有一定的反射作用，同时扩散的核素粒子受沉降等作用，核扩散物质又不可能被全部反射回去，因而核辐射物质只能是部分反射回大气，为了便于描述和模拟，设地面反射系数为。这样进入大气的核辐射物质可以看成是两个部分：一是从泄漏源直接扩散到空间A点；二是从地面反射进入空间A点（见图3）图3 连续点源扩散地面部分示意图从几何物理学分析可知，通过地面反射进入大气的核扩散物质相当于虚拟泄漏源泄露的核辐射物质在原来空间的一个浓度叠加。如果所设核泄漏点源在距有效地面的高度为的地方，本文求空间任一点的浓度值，则实际泄漏源对点的影响部分可用来表示。在考虑反射系数后，虚泄漏源对点的影响部分可用来表示，于是（10）式所得模型可以修正为： （16）5.2问题二模型的建立与求解5.2.1.1模型二的建立问题二，针对在有风速时，选取高斯连续点源扩散模型进行分析，这也是在所有有风速时，放射性物质扩散模型中最常见也最方便的一种泄漏方式。考放射性气体云团在大气中实际迁移和扩散的数值计算基本上可分为二步：第一步，根据大气动力学理论进行所关注区域中风场的计算，其理论基础是大气运动方程、连续性方程、状态方程、热力学方程和水汽方程构成的基本方程组；第二步，进行已知风场中放射性气体云团迁移和扩散的计算，可采用类似于处理大气污染的方法，假设放射性气体云团不影响大气流体速度和温度，求解出放射性气体云团的连续性方程。高斯连续点源扩散模型有边界点源，这里设泄漏源有效高度为，取其地面投影为坐标原点，轴指向风向，继续考虑地面反射作用，图4 边界点源示意图可得到高斯连续点源泄漏的浓度分布7： （17）其中，为源强，；为污染物质量浓度，；分别为用浓度标准差表示的轴上的扩散参数；为泄漏高度的平均风速，m/ s；为泄漏有效高度，。考虑到动力抬升的作用，将有效高度（即有效烟云高度）。图5 动力抬升示意图因此（17）式可以写成为： （18）5.2.1.2对高斯浓度计算模型的修正利用标准的高斯模型，计算大气中放射性核素云团的扩散，目前已有较多的研究，而实际的核素扩散过程还存在着粒子的重力沉降、雨洗沉积以及核素衰变等的因素对浓度分布的影响。因此，有必要对高斯浓度扩散模式进行修正，方可较为准确、真实地反映实际的核素扩散规律。由于地面风速对大气扩散的影响起着至关重要的作用，不同的风速其浓度计算方法有较大差异，因此，下面的分析从有风（大于 1. 5 m/s）、小风（小于1. 5m/s大于0. 5m/s）和静风（小于0. 5m/ s），三种情况分别论述。a：考虑干沉积时的连续点源扩散粒径大于10 m 的粒子有明显的重力沉降，粒子的沉降速度取决于空气阻力和重力平衡，可用斯托克斯公式表示： （19）式中，：粒子密度， ；：重力加速度，9. 806 5；：粒子直径，;：空气的动力粘性系数, 可取；：沉降速度，含碘放射性核素的干沉积速度为= 1. 1 cm/ s, 该值也可由式( ) 计算得到。因为在扩散过程中同时有重力沉降的位移迭加到羽流中心线上，中心线就会向下倾斜， 所有粒子相当于在下倾的中心线上扩散。该类扩散和沉降的迭加可认为是羽流运行过程中，实源以的速度向下移动，在处向下移动的高度为，即源高由降到了。实际上，由于大气湍动及其他动力作用， 地面不是全吸收表面，应考虑地面的反射作用。但毕竟存在粒子的沉降作用，又不是全反射，因此，需对反射项乘以一个反射系数，由于反射项的有效源高度也变成了，故相应的浓度计算公式为：（20）式( 20) 可得到经过干沉积作用后的浓度分布结果。式( 20) 中, 反射系数需外部给定， 通常对于放射性核素可取0.5。b: 小风（小于1. 5m/s大于0. 5m/s）和静风（小于0. 5m/ s）时的计算参考小风和静风时的浓度计算方法，再利用部分反射倾斜烟云模式，可得出该情况下的浓度分布：（21）通常按照每小时排放6个烟团进行积分，即可达到需要的精度。c:地面的干沉积量根据扩散理论和动量传递的普朗特理论，可以导出干沉积的地面沉积量为： （22）式中，前面各种情形计算的地面污染物的浓度；干沉积速度，地面干沉积率，。d: 雨洗作用降雨对烟羽中的颗粒物及气溶胶具有清洗作用，可溶性气体与蒸汽亦可溶于雨水中，降雨过程造成的这类湿沉积是导致放射性气溶胶和气体向地面沉积的另一重要机制。通常以冲洗系数，描述降雨对烟羽中污染物清洗作用的大小。与雨强的关系可以表达为： （23）式中，为雨强(mm/ h) ；a，b为经验系数。式中，按释放物质为含碘、不含碘情况分别取值。对含碘物质，取；对于不含碘物质，取。对于湿沉积导致的烟羽耗减，可采用湿沉积耗减因子对源强进行修正，有 （24）e：放射性衰变的影响除了干沉积和湿沉积外，放射性物质的衰变也是影响大气中核素浓度分布的主要原因。由于放射性物质服从简单的衰变规律，其浓度随时间的变化可由下式计算： （25）式中， 为初始浓度；衰变常数；经过的时间；可由无衰变的浓度公式计算。参照前面的模式，同样采用衰变耗减因子来对源项进行修正。经过理论推导，可以得出： （26）其中，放射性核素的半衰期。如以上情况都存在，应同时考虑源项修正，即式(24) 和式( 26) 相乘后代替5.2.2模型二的求解5.2.2.1模型参数的确定在事故过程模型中所需参数的选取及确定十分重要，通常情况下气象参数的选取是利用该地区多年气象资料，采取工业安全与环保统计的方法进行有关参数的确定，而其他扩散参数是以实际测定为基础的。a:大气稳定度的计算根据国家标准(GB/ T 13201 - 1991) )制定地方大气污染物排放标准的技术方法的规定，划分大气稳定度的级别，共分为6级A F，A为极不稳定；F为极稳定。首先，根据释放源所在地的经度和纬度以及泄漏的日期和时间计算当时的太阳高度角h0；然后，由太阳高度角h0和云量查出太阳辐射等级；最后，再根据地面风速确定当时的大气稳定度，计算细节可参考文献。表 1 大气稳定度的级别参考表地面风速白天太阳辐射阴天的白天或夜间有云的夜晚强中弱薄云遮天或低云0.5云量0.46CDDDDDb：扩散参数的计算有风时的扩散参数，的确定采用Briggs给出一套扩散参数，如表1和表2所示。表2 Briggs 扩散参数（开阔平原田野）大气稳定度ABCDEF表3 Briggs 扩散参数（城市）大气稳定度A-BCDE-F5.2.2.2模型的求解图6 模型二计算流程图5.3问题三模型的建立与求解5.3.1模型三的建立问题三，针对在有风速时，选取高斯连续点源扩散模型进行分析在给出的上风和下风公里的距离处在结合，在模型二的基础上考虑浓度等效圆的情况。在模型二里面将式（18）中的风速在上风和在下风不同情况下与传播速度之间的比较的分析。a:在上风的情况下，此时满足（）此时若时，即风速小于自然扩散速度，此时的放射性物质是无法到达上风口的，则必须满足风速大于自然扩散速度，仍然根据模型二可以求解，只是将式（18）中的用来代入。b:在上风的情况下，此时无较大的影响和要求此时不论风速小于自然扩散速度，还是风速大于自然扩散速度，此时的放射性物质都可到达下风口的，则仍然根据模型二可以求解，步骤同5.2.2.图7 核泄漏对上下风在距离上的影响5.4问题四的求解日本于2011年3月15日下午，世界气象组织和国际原子能机构北京区域环境紧急响应中心分析认为：日本中北部区域在中低层大气中的风向由西南风转为西北风；高空大气主要以偏西风气流为主，近期由于降水发生，有利于核物质沉降，影响范围缩小。未来三天（16日至18日），日本核电站核泄漏产生的放射性污染物主要影响区域为日本中部、北部及其以东的北太平洋区域。在结合模型二、模型三的条件下，在参阅整理大量气象、地理、新闻资料，选择我国东海岸山东半岛作为研究对象，综合考虑对应海域平均风速及风向、地理距离、海水对放射性物质扩散的部分反射系数等因素，并通过C+编程模拟计算，预测出放射性核物质将经过6.5天到达我国东海岸山东半岛，且131I浓度值为：0.100与，与实际情况比较吻合。6结果分析与检验由于在前四问题中均无具体的数据，所以对于结果，在检验上存在一些麻烦，但是在模型三预测上下风的浓度具有较好的实际意义，而在模型四中，对于预测出浓度的变化和对我国的影响也较为合理，故该模型具有可靠性和参考性，可以进行推广和其它方面的应用。7模型的优缺点与改进方向优点：该方法可计算核事故中连续点源和瞬时点源在不同气象、地形条件下的浓度分布。该结果在核事故的应急救援过程中，对救援人员划定警戒区和确定周围居民的疏散范围具有重要意义，并可为制定救援方案和应急决策提供科学依据。放射性云团在空中迁移和扩散提供浓度的定量描述，对可能发生的核事件的放射性核素浓度监测及监测时间范围提供相关信息。缺点：模型在建立时只考虑了两种下垫面情况，即农田和城市，实际的扩散情况可能因地形条件的复杂化而有所改变。由于扩散是一个非常复杂、影响因素众多的过程，上述方法还有一定的局限性，有待于进一步检验。另外，该方法只给出了浓度的计算结果，没有换算成人的吸收剂量，今后将对其计算过程和参数的选取进行多方面验证，使其进一步完善。参考文献1金勇进,统计学,北京：中国人民出版社,2010。2梅长林,周家良,实用统计方法,北京:科学出版社,2002。3罗艾民,魏利军.有毒重气泄露安全距离数值方法J.中国安全科学学报, 2005,15（8）: 98-100 4周燕等,Matlab在统计与工程数据分析中的应用,北京:电子工业出版社,2010。5丁信伟,王淑兰,徐国庆.可燃及毒性气体泄漏扩散研究综述.化学工业与工程,1999,2(16):118122.6胡世明.气体释放源的三维瞬态重气扩散研究.劳动保护科学技术,2O02,3(20):28307蒋军成,潘旭海.描述重气扩散的一种新模型.南京工业大学学报,2OO2,I(24):41463附 录附录一：模型二在上下风L公里的浓度求解源代码#include#include#includeusing namespace std;#define e 2.718281828459 #define PI 3.1415926void qiuxy(char wending ,double x,double &dy,double &dz) switch(wending) case A: dy=0.22*x/(sqrt(1+0.0001*x), dz=0.2*x; break; case B: dy=0.16*x/(sqrt(1+0.0001*x), dz=0.12*x; break; case C: dy=0.11*x/(sqrt(1+0.0001*x), dz=0.08*x/(sqrt(1+0.0002*x); break; case D: dy=0.08*x/(sqrt(1+0.0001*x), dz=0.06*x/(sqrt(1+0.0015*x); break; case E: dy=0.06*x/(sqrt(1+0.0001*x), dz=0.03*x/(1+0.0003*x); break; case F: dy=0.04*x/(sqrt(1+0.0001*x), dz=0.016*x/(1+0.0003*x); break;void main() double c1,c2,x,z,t,p0,H,dx,dy,dz,a,s,k,I,Ts; char daqi;cout请按以下格式输入：endl;cout下风口L处坐标（米）：x;cout地点的纵坐标（米）：endl; cout0 endl; cout地点的铅直坐标（米）：z;cout扩散时间（秒）：t;cout地面反射系数：a; cout放射性气体扩散速度（m/s）：s;cout风速（m/s）：k; cout雨水强度（毫米/小时）：I;cout核物质衰变周期（天）：Ts; cout初始浓度（千克/立方米）：p0;cout有效高度（米）：H;cout大气稳定度等级：daqi;qiuxy(daqi,x,dy,dz);dx=dz; c1=p0/(2*PI*(s-k)*dy*dz