空气污染数学建模

1A.污染气体的传播扩散摘要钢铁生产排放的污染气体是造成雾霾的重要原因之一，研究污染气体的扩散特征，正确模拟污染气体的扩散过程，能够为钢铁生产集团提出更好的治理管理措施，具有实际意义。针对问题一：污染气体的排放速度为300m/s,在不考虑风向风速及高度影响的情况下，此问题即为二维平面的连续点源扩散问题，由此在二维xoy平面上建立连续点源扩散方程模型()(,,)txxyyuuufxyt，其中为气体扩散系数，本文中取为常数10，f(x,y,t)为污染气体的排放速度，在本文中恒为300m/s；对上述偏微分方程模型，本文采用ADI法（Alternatingdirectionimplicit，交替方向隐式法）求解出迭代格式，利用MATLAB编程，求出模型一的数值解，并得到任意时刻污染气体的浓度分布情况。通过SPSS软件，对附件一所给的原始实际数据与模型一求解得到的模拟值进行显著性检验，检验结果显示该模型与实际情况吻合。针对问题二：考虑风向风速对污染气体扩散过程的影响时，在基于对问题一求解的基础上，在模型一的扩散方程模型中加入风向风速的平流项，由此得到有风情况下的模型),,()(21tyxfuuuuuyyxxyxt，其中12，分别为风速在x,y方向的分量；对此模型同样采用ADI法求出迭代格式，利用MATLAB编程，求出模型二的数值解，并得到任意时刻污染气体的浓度分布情况。通过SPSS软件，对附件二所给的原始实际数据与模型二求解得到的模拟值进行显著性检验，检验结果显示该模型与实际情况吻合。针对问题三：考虑有风时增加高度的影响，此问题即为三维空间的污染气体扩散问题，考虑到三维模型的编程复杂度，而且污染气体的扩散在xoy平面上各向同性，可以将污染气体在y方向的扩散等价为在x方向上的扩散，此时便只需要建立xoz平面上的扩散模型。在基于对问题二求解的基础上，在模型二的扩散方程中增加高度项，由此得到模型三为12()2,,txxxzzuzuuufxzt，其中为z方向的扩散系数；对该扩散方程同样采用ADI法求出迭代格式，利用MATLAB编程，求出模型二的数值解，并得到任意时刻污染气体的浓度分布情况。关键词：污染气体扩散方程ADI法数值解2一、问题重述目前，治理雾霾是人们最为关心的热点问题之一。中国社科院发布的《气候变化绿皮书》中提及，雾霾形成的原因里，重工业、车辆尾气、土方施工都榜上有名，其中钢铁生产也是造成雾霾的重要原因之一。某钢铁生产集团烟囱污染气体的排放对周边地区大气污染的影响非常大，为了提出更好的治理管理措施，需要对其污染气体扩散的特征进行分析。现在，我们需要在三种情况下考虑污染气体的扩散过程：1.在不考虑风向和高度影响的情况下，建立模型，模拟某钢铁生产集团的烟囱排放污染气体的扩散过程，假设烟囱的排放速度为300m/s。2.考虑风向为东北风，平均风速0.6m/s的情况下，模拟污染气体的传播扩散过程。3.在考虑风向的基础上增加高度的影响，建立模型，模拟污染气体的传播扩散过程。4.基于上述模型结论，给该钢铁生产集团提供一个污染气体治理建议报告。二、问题分析钢铁生产集团烟囱污染气体的排放对周边地区大气污染的影响非常大，为了提出更好的治理管理措施，需要对其污染气体扩散的特征进行分析。污染气体扩散是与时间、空间相关的连续性问题，本文利用微分方程建立模型，通过ADI法（Alternatingdirectionimplicit，交替方向隐式法）求解微分方程的迭代式，利用MATLAB编程求其数值解来模拟扩散过程。对于问题一，在不考虑风速和高度影响的情况下，污染气体的扩散问题即为二维平面点源的扩散问题；考虑到烟囱的排放速度为300m/s，需要在二维扩散方程中加入有源项，由此建立扩散模型。为检验模型的准确性，需利用SPSS软件，将扩散方程的数值解与实际数据进行显著性分析。对于问题二，在问题一的基础上，考虑了风向风速的影响，需要在模型一中加入风向风速对应的平流项，由此建立模型二。为检验模型的准确性，同样需要利用SPSS软件，将模型二的数值解与实际数据进行显著性分析。对于问题三，基于对问题一、二的建模思路，首先建立不考虑风速风向只考虑高度时的模型，此时应分析扩散系数与高度的关系，建立模型三；然后在模型三的基础上考虑风速风向的影响，此时应分析风速与高度的关系，建立模型四。为使模型更能符合实际情况下污染气体的扩散过程，本问题还考虑了大气静力稳定度对扩散过程的影响。综上所述，本问题可以看成是求解偏微分方程中扩散方程的数值解问题。三、模型假设与符号说明1.假设烟囱的排放速度恒为300m/s；2.假设污染气体中只含有同一类污染物，污染气体的扩散系数相同；3.假设不考虑气体受到的重力和浮力；4.假设在整个扩散过程中污染气体不发生沉降、分解，不发生化学反应；5.假设地面以及地标地物对气体无吸收；36.假设在有风情况下，风向水平，风速风向恒定；u:污染物浓度f:污染物排放速度v:风速:水平方向扩散系数1：风速在x方向的分量2：风速在y方向的分量：垂直方向扩散系数四、对问题一的分析与建模4.1问题分析问题一在不考虑风速和高度影响下，烟囱以恒定速度排放污染气体，在无限大平面上，将烟囱口看成是点源，此问题即为二维平面点源的扩散问题，建立二维扩散方程模型，通过ADI法求其数值解，即可得到任意时刻污染气体的浓度分布。4.2模型一的建立本文将烟囱口当作二维平面上的点源来考虑，分析附件allresults_1所给的数据，可以得出的结论有：污染气体扩散的平面图为1000×1000（单位：m）的正方形区域；以点（500,500）为中心，半径为1m的领域内污染气体的浓度最高，据此确定烟囱口的位置坐标为（500,500）。假设污染源为连续点源，建立有源项的二维扩散方程：1000,01000,0,),,()(yxtyxfuuuyyxxt（1）其中，为扩散系数，t为时间，f为污染源排放污染气体的速度，tu为浓度对时间的偏导，xxu为浓度对x的二阶偏导，yuy为浓度对y的二阶偏导。现在考虑初边值条件，初始时刻，平面上各点处浓度均为0，即有(,,0)0,,[0,1000]uxyxy边界条件取延拓边界条件，即有(0,,)(1000,,)(,0,)(,1000,)0,,[0,1000]uytuytuxtuxtxy根据题意知，在污染源处，污染物排放速度为300m/s，假设烟囱口是以点（500,500）为中心，1m为半径的圆形区域，在1000×1000的平面内，上述圆形区域任然可以当作点源来处理，在圆形区域内有源项为300，在圆形区域外有源项为0，即有422300,(-500)(-500)1(,,)0,,xyfxytxy其他其他由此得到在不考虑风向和高度影响时连续源的扩散模型()(,,),0,1000(0,,)(1000,,)0,,[0,1000](,0,)(,1000,)0,,[0,1000](,,0)0,,[0,1000]txxyyuuufxytxyuytuytxyuxtuxtxyuxyxy（2）另外，根据《中国大百科全书》大气扩散卷，选定扩散系数sm/102。4.3模型一求解本文采用ADI法求解（1）式数值解，将（1）式转化为ADI格式为111122221,,1,,1,,1,2222/2nnnnnnnnjkjkjkjkjkjkjkjkjkuuuuuuuuaaftxy（a）1111122221,,1,,1,,1,2222/2nnnnnnnnjkjkjkjkjkjkjkjkjkuuuuuuuuaaftxy（b）对（a）式进行转换，令：22,ytarxtaryx并将yxrr,带入（a）式中，得：1112221,,1,,1,,1,11()(1)()22111()(1)()222nnnxjkxjkxjknnnyjkyjkyjkjkrurururururutf(3)将(3)式写成矩阵方程组形式：51111221112211122111221xxxxxxxxxxxxrrrrrrrrrrrr12,112,212,312,212,112,njnjnjnjKnjKnjKuuuuuuyyyy1111221112211122111221yyyyyyyyrrrrrrrrrrrr,1,2,3,2,1,njnjnjnjKnjKnjKuuuuuuJ,...,3,2,1j同样对（b）式进行转换，令22,ytarxtaryx，并将yxrr,带入（a）式中，得：kjnkjxnkjxnkjxnkjynkjynkjytfurururururur,21,121,21,111,1,11,21)21()1()21()21()1()21((4)将(4)式写成矩阵方程组形式：61111221112211122111221xxxxxxxxxxxxrrrrrrrrrrrr121,122,123,122,121,12,nknknknJknJknJkuuuuuu1111-22111-22111-22111-221xxxxxxxxxxxxrrrrrrrrrrrr1,2,3,2,1,,nknknknJknJknJkuuuuuuKk3,2,1根据（3）（4）的迭代格式及（2）中初边值条件，运用MATLAB编程，可以得到任意时刻的污染气体浓度分布，其中截取60s、600s、1800s和3600s时浓度分布如图1.图1（a）60s时浓度分布图图1（b）600s时浓度分布图7图1（c）1800s时浓度分布图图1（d）3600s时浓度分布图1不同时刻浓度分布图从图1不同时刻浓度分布图可以看出：开始时（t=60s），烟囱口处污染气体的浓度近似为800个单位，远远高于周围其他区域，整个浓度分布区域近似为一条垂直于xoy平面的直线；随着时间的推移，到t=600s时，烟囱口处污染气体的浓度近似为1100个单位，且周围其他区域内的浓度有所增大，整个浓度分布区域由60s的直线状变为呈上尖下圆的塔状；当t=1800s时，烟囱口处污染气体的浓度近似为1300个单位；当t=3600s时，烟囱口处污染气体的浓度已达到1700个单位左右，从“塔”的形状来看，“塔尖”在不断地增高，“塔底”在不断地向外扩展，到3600s时已扩散到半径为800米的圆形区域内。4.4模型一的检验首先通过MATLAB画图，比较模型一的浓度分布图与用实际数据画出的浓度分布图的区别，截取120s和3600s的浓度对比图如下：图2（a）120s时模拟值与实际数据浓度分布对比图图2（b）3600s时模拟值与实际数据浓度分布对比图图2不同时刻模拟值与实际数据浓度分布对比图由图2（a）（b）两幅图可以看出：由模型一得到的浓度分布图与用实际数据画出的8浓度分布图在浓度大小和浓度分布形状上近似一致。为了更进一步的分析两者之间的误差，现在运用SPSS进一步做模型检验。取120s和3600s时由模型一得到的浓度分布数据与实际数据进行显著性检验。首先画出模拟值与实际数据的比较图如下：图3（a）120s模拟值与实际数据比较图图3（b）360