污染气体扩散文档.docx
《污染气体扩散文档.docx》由会员分享，可在线阅读，更多相关《污染气体扩散文档.docx（35页珍藏版）》请在冰豆网上搜索。
污染气体扩散文档
污染气体扩散模拟程序设计
目录
1引言1
2大气扩散模式2
2.1大气扩散模式概述2
2.2高斯大气扩散模式2
2.3高斯计算公式4
2.3.1连续泄露高斯模型4
2.3.2瞬时泄漏高斯模型4
2.4高斯扩散方程的求解4
2.4.1大气稳定度的确定4
2.4.2高架源回归系数（Py、qy、Pz、qz）确定5
2.4.3设计流程图6
3程序实现步骤6
3.1部署界面6
3.2定义及初始化相关参数8
3.2.1添加相关命名空间8
3.2.2定义相关参数8
3.2.3初始化相关参数9
3.2.4获取采样范围13
3.2.5自定义Polygon类16
3.3单源瞬时大气污染的高斯扩散模拟18
3.3.1添加Timer控件18
3.3.2开始模拟19
3.3.3在内存中创建featureLayer22
3.3.4ContrainAngle函数25
3.3.5瞬时泄露高斯计算函数25
3.3.6反距离加权内插26
3.3.7分级渲染27
4总结28
1引言
有毒或易燃易爆危险物质,一旦发生泄漏事故,将对周围环境和人员造成巨大危害。
大量的易挥发、易燃性物质迅速释放到大气中,会形成气云并逐渐扩散。
对于易燃易爆物质假如气云在尚未稀释到低于燃烧下限前被引燃,就可能发生气云爆炸或闪火,并引起大范围的破坏对于有毒物质,在扩散范围内蒸气云团的浓度稀释至安全浓度前能引起中毒。
本文介绍如何在.net平台下，利用ArcGISEngine完成基于高斯模型的瞬时污染气体扩散软件。
2大气扩散模式
2.1大气扩散模式概述
大气扩散模式是对大气污染物散布过程的模拟。
按污染源的几何特征可分为电源、线源、面源扩散模式；按污染源排放特性分有正常工况模式和非正常排放模式；按假设条件不同分有烟流扩散模式和烟团扩散模式；按模式适用范围不同分有短距离扩散模式和中长距离输送模式等。
具体操作时应根据当地地形、气候等环境要素以及污染源亲何况正确选取并适当修正。
主要的大气扩散模式有：
（1）高架电源扩散模式（静风模式、有效源高等）
（2）面源扩散模式（虚电源模式、窄烟云模式：
在多面源情形下计算效果较好）
（3）线源扩散模式（电源求和法、风向与线源平行时扩散模式）
2.2高斯大气扩散模式
高斯大气扩散模式属于高架电源的烟羽大气扩散模式，它是计算释入大气中的气载污染物下风向浓度的应用最广的方法。
高斯扩散模式所描述的扩散过程主要有以下几种：
下垫面平坦、开阔、性质均匀,平均流场平直、稳定，不考虑风场的切变；扩散过程中,污染物本身是被动、保守的,即污染物和空气无相对运动,且扩散过程中污染物无损失、无转化,污染物在地面被反射；扩散在同一温度层集中发生,平均风速>1.0m/s;适用范围一般<10～20km。
在均匀、定常的湍流大气中污染物浓度满足正态分布，但实际大气一般不满足均匀、定常条件，因此高斯模式应用于下垫面均匀平坦、气流稳定的小尺度扩散问题更为有效。
高斯模式的有关假定：
坐标系
原点为高架点源排放点在地面的投影点，x轴正方向为平均风向，y为横风
向，在水平面上垂直于x轴，z轴垂直于水平面xoy，向上为正向。
四点假设
a．污染物浓度在y、z风向上分布为正态分布
b．全部高度风速均匀稳定
c．源强是连续均匀稳定的
d．扩散中污染物是守恒的（不考虑转化）
高斯扩散模式的坐标系：
高斯烟流的浓度分布：
高斯烟流中心线上的浓度分布
2.3高斯计算公式
2.3.1连续泄露高斯模型
对于在恒定气象条件（指风向、风俗、大气稳定度不随时间而变）下的高架电源的连续排放，在考虑了烟羽在地面的全反射后，下风向任一点的污染物浓度c（x,y,z,H）可由下式计算：
}C（x,y,z）：
下风向某点（x,y,z）处的空气中污染物的浓度（mg/m3）；
}x：
下风向距离（m）；
}y：
横截风向距离（m）；
}z：
距地面的高度（m）；
}q：
气载污染物源强，即释放率（mg/s）；
}u：
排气筒出口处的平均风速（m/s）；
}H：
排气筒距地面的几何高度（m）；
}δx,δy,δz：
扩散系数
2.3.2瞬时泄漏高斯模型
其中t为时间，单位为秒，其余参数与连续泄露高斯模型中的含义一样。
注:
该程序使用瞬时泄露高斯模型。
2.4高斯扩散方程的求解
2.4.1大气稳定度的确定
大气稳定度的确定根据五类地面风速、三类日间日射和两类夜间云量把扩散天气分为6类，即强不稳定、不稳定、弱不稳定、中性、较稳定和稳定。
分别用英文字母A、B、C、D、E和F表示。
在国标“制订地方大气污染排放标准的技术原则和方法”（GB3840-83）与“环境影响评
价技术原则”（HJ/T2.1-93）中，建议采用下属修订的帕斯奎尔稳定度分类方法。
首先由云量与太阳高度角按表1查出太阳辐射等级数，再由太阳辐射等级数与地面风速按表2查找稳定度等级。
表1太阳辐射等级数
表2大气稳定度等级
2.4.2高架源回归系数（Py、qy、Pz、qz）确定
国际原子能机构（IAEA）根据相关实验数据，推荐了一组适用于大粗糙度（Z0≈1m）和高架源的弥散系数公式，如表3，据此求得相应地扩散参数。
表3与排放高度有关的扩散参数
2.4.3设计流程图
其中1回归系数由上述三个表可得；2采样点位置，根据风向与污染源位置计算出可能受影响范围，然后设置采样间隔进行采样；3污染源浓度默认为500,；4高斯模型即高斯计算公式。
3程序实现步骤
3.1部署界面
在VisualStudio2005或其它较高版本下，“新建项目”－“ArcGIS”－“Engine”，选择“MapControlApplication”,按照下图，部署界面。
其中右边面板：
将“播放速度”处的trackBar控件的属性设置为：
Maximum：
5，Minimum:
1。
3.2定义及初始化相关参数
3.2.1添加相关命名空间
usingESRI.ArcGIS.esriSystem;
usingESRI.ArcGIS.Carto;
usingESRI.ArcGIS.Controls;
usingESRI.ArcGIS.ADF;
usingESRI.ArcGIS.SystemUI;
usingESRI.ArcGIS.Geometry;
usingESRI.ArcGIS.Display;
usingESRI.ArcGIS.Geodatabase;
usingESRI.ArcGIS.Geoprocessing;
usingESRI.ArcGIS.SpatialAnalyst;
usingESRI.ArcGIS.GeoAnalyst;
usingESRI.ArcGIS.DataSourcesRaster;
usingESRI.ArcGIS.DataSourcesGDB;
3.2.2定义相关参数
publicstructDiffusedParameter//自定义一个大气扩散参数结构体
{
publicdoublePy;
publicdoublePz;
publicdoubleQy;
publicdoubleQz;
}
//初始化扩散参数表，此时的0、1、2、3、4、5对应大气稳定度等级A、B、C、D、E、F
publicDiffusedParameter[]diffusedParameter50;
publicDiffusedParameter[]diffusedParameter100;
publicDiffusedParameter[]diffusedParameter180;
constdoublePI=3.1415926;
privateintradius=-1;//辐射等级
privateintatmosphereStability=-1;//初始化大气稳定度等级
privatedoublePy,Pz,Qy,Qz;//计算中使用的扩散参数
privatedoubleQ=500;//气载污染物源强，即释放率（mg/s）
IEnvelopeenvelope=newEnvelopeClass（）;//初始5公里范围
IPointwuranyuan=newPointClass（）;//定义污染源点
doubleinitialscale=1;//初始比例
boolisDoubleClick=false;//定义是否双击MapControl的标示变量
IPointCollectionpPointCollection;
IPolygonpPolygon;
IEnvelopepEnvelope;
3.2.3初始化相关参数
3.2.3.1选中右边面板所有的RadioButton按钮，在radioButton1中响应其Click事件，即可以实现对所有RadioButton控件的Click事件的响应
privatevoidradioButton1_MouseClick（objectsender,MouseEventArgse）
{
RadioButtonradioButton=（RadioButton）sender;
radius=Convert.ToInt32（radioButton.Text）;//将radioButton的Text值赋给radius变量
}
3.2.3.2获取污染源位置
在axMapControl1的OnDoubleClick添加如下代码，以获取污染源的位置。
IGraphicsContainerpGra;
privatevoidaxMapControl1_OnDoubleClick（objectsender,IMapControlEvents2_OnDoubleClickEvente）
{
if（this.axMapControl1.LayerCount==0）{MessageBox.Show（"请先加载地图！
","提示",MessageBoxButtons.OK,MessageBoxIcon.Warning）;return;}
isDoubleClick=true;
wuranyuan.X=e.mapX;
wuranyuan.Y=e.mapY;
IMarkerElementmarkerelement=newMarkerElementClass（）;
IElementpElement=markerelementasIElement;
pElement.Geometry=wuranyuan;
pGra=this.axMapControl1.MapasIGraphicsContainer;
pGra.DeleteAllElements（）;
pGra.AddElement（pElement,1）;
this.axMapControl1.ActiveView.PartialRefresh（esriViewDrawPhase.esriViewGraphics,null,null）;
}
3.2.3.3添加initialvar（）函数，用于初始化相关参数
privatevoidinitialvar（）
{
pPointCollection=null;
pPolygon=null;
pEnvelope=null;
//根据地面风速和太阳辐射等级计算大气稳定度
#region
doublespeed;
try
{
speed=Convert.ToDouble（this.windspeed.Value）;
}
catch
{
MessageBox.Show（"风速填写错误，请查看！
"）;
return;
}
if（speed<0）
{
MessageBox.Show（"风速填写错误，请查看！
"）;
return;
}
switch（radius）
{
case-2:
if（speed<=2.9）
atmosphereStability=6;//此时的atmosphereStability的值1、2、3、4、5、6对应大气稳定度等级A、B、C、D、E、F
elseif（speed>2.9&&speed<=4.9）
atmosphereStability=5;
elseif（speed>4.9）
atmosphereStability=4;
break;
case-1:
if（speed<=2.9）
atmosphereStability=5;
elseif（speed>2.9）
atmosphereStability=4;
break;
case0:
atmosphereStability=4;
break;
case1:
if（speed<=1.9）
atmosphereStability=2;
elseif（speed>1.9&&speed<=4.9）
atmosphereStability=3;
elseif（speed>4.9）
atmosphereStability=4;
break;
case2:
if（speed<=1.9）
atmosphereStability=1;
elseif（speed>1.9&&speed<=4.9）
atmosphereStability=2;
elseif（speed>4.9&&speed<=5.9）
atmosphereStability=3;
elseif（speed>5.9）
atmosphereStability=4;
break;
case3:
if（speed<=2.9）
atmosphereStability=1;
elseif（speed>2.9&&speed<=4.9）
atmosphereStability=2;
elseif（speed>4.9&&speed<=5.9）
atmosphereStability=3;
elseif（speed>5.9）
atmosphereStability=4;
break;
}
#endregion
//根据大气稳定度和释放高度计算回归系数Py,Pz,Qy,Qz
#region
if（height.Text==""）
{
MessageBox.Show（"请设置释放高度参数！
"）;
return;
}
if（height.Text.Trim（）=="50"）
{
Py=diffusedParameter50[atmosphereStability-1].Py;
Pz=diffusedParameter50[atmosphereStability-1].Pz;
Qy=diffusedParameter50[atmosphereStability-1].Qy;
Qz=diffusedParameter50[atmosphereStability-1].Qz;
}
if（height.Text.Trim（）=="100"）
{
Py=diffusedParameter100[atmosphereStability-1].Py;
Pz=diffusedParameter100[atmosphereStability-1].Pz;
Qy=diffusedParameter100[atmosphereStability-1].Qy;
Qz=diffusedParameter100[atmosphereStability-1].Qz;
}
if（height.Text.Trim（）=="180"）
{
Py=diffusedParameter180[atmosphereStability-1].Py;
Pz=diffusedParameter180[atmosphereStability-1].Pz;
Qy=diffusedParameter180[atmosphereStability-1].Qy;
Qz=diffusedParameter180[atmosphereStability-1].Qz;
}
#endregion
//根据污染源位置与风向确定可能受污染范围
//初始比例
initialscale=10;
IEnvelopemapenvelope=this.axMapControl1.Extent;
if（（wuranyuan.X+5000/initialscale）>mapenvelope.XMax）{envelope.XMax=mapenvelope.XMax;}
else{envelope.XMax=wuranyuan.X+5000/initialscale;}
if（（wuranyuan.X-5000/initialscale）else{envelope.XMin=wuranyuan.X-5000/initialscale;}
if（（wuranyuan.Y+5000/initialscale）>mapenvelope.YMax）{envelope.YMax=mapenvelope.YMax;}
else{envelope.YMax=wuranyuan.Y+5000/initialscale;}
if（（wuranyuan.Y-5000/initialscale）else{envelope.YMin=wuranyuan.Y-5000/initialscale;}
//采样点图层
pEnvelope=newEnvelopeClass（）;
pPolygon=GetInterpolateArea（）;//得到采样范围
pPointCollection=（IPointCollection）pPolygon;
}
3.2.4获取采样范围
publicIPolygonGetInterpolateArea（）
{
//定义变量
#region
IPointpt=newPointClass（）;
pt=wuranyuan;
//if（winddirection.Text==""）{MessageBox.Show（"请输入风向!
"）;returnnull;}
intdirection=Convert.ToInt32（winddirection.Text）;//风向
doubleangle=（direction+90）%180;//与风向垂直方向的夹角
doubletan=Math.Tan（PI*angle/180）;//与风向垂直的斜率
IPointpt1=newPointClass（）;
intpt1Flag=0; //记录第一个交点在那条边上1234分别表示上右下左边
IPointpt2=newPointClass（）;
intpt2Flag=0;//记录第二个交点在那条边上
doubletempX,tempY;
#endregion
//依次判断,上边
#region
tempY=envelope.YMax;
tempX=（tempY-pt.Y）/tan+pt.X;
if（tempX>=envelope.XMin&&tempX<=envelope.XMax）
{
pt1.X=tempX;
pt1.Y=tempY;
pt1Flag=1;
}
#endregion
//右边
#region
tempX=envelope.XMax;
tempY=tan*（tempX-pt.X）+pt.Y;
if（tempY>=envelope.YMin&&tempY<=envelope.YMax）
{
if（pt1Flag==0）
{
pt1.X=tempX;
pt1.Y=tempY;
pt1Flag=2;//当该点为第一个交点时，即标识第一个交点在第2条边上
}
else
{
pt2.X=tempX;
pt2.Y=tempY;
pt2Flag=2;//若该交点为第二个交点，则标识第二个交点在第2条边上
}
}
#endregion
//下边
#region
tempY=envelope.YMin;
tempX=（tempY-pt.Y）/tan+pt.X;
if（tempX>=envelope.XMin&&tempX<=envelope.XMax）
{
if（pt1Flag==0）
{
pt1.X=tempX;
pt1.Y=tempY;
pt1Flag=3;
}
elseif（pt2Flag==0）
{
pt2.X=tempX;
pt2.Y=tempY;
pt2Flag=3;
}
}
#endregion
//左边
#region
tempX=envelope.XMin;
tempY=tan*（tempX-pt.X）+pt.Y;
if（tempY>=envelope.YMin&&tempY<=envelope.YMax）
{
if（pt1Flag==0）
{
pt1.X=tempX;
pt1.Y=tempY;
pt1Flag=4;
}
else
{
pt2.X=tempX;
pt2.Y=tempY;
pt2Flag=4;
升级会员 