A.污染气体的传播扩散摘要钢铁生产排放的污染气体是造成雾霾的重要原因之一，研究污染气体的扩散特征，正确模拟污染气体的扩散过程，能够为钢铁生产集团提出更好的治理管理措施，具有实际意义。

针对问题一：污染气体的排放速度为300m/s,在不考虑风向风速及高度影响的情况下，此问题即为二维平面的连续点源扩散问题，由此在二维xoy 平面上建立连续点源扩散方程模型，其中为气体扩散系数，本文中取为常数10，f(x,y,t ) 为污染气体的排放速度，在本文中恒为300m/s ；对上述偏微分方程模型，本文采用ADI法（Alternating direction implicit，交替方向隐式法）求解出迭代格式，利用MATLAB编程，求出模型一的数值解，并得到任意时刻污染气体的浓度分布情况。

通过SPSS软件，对附件一所给的原始实际数据与模型一求解得到的模拟值进行显着性检验，检验结果显示该模型与实际情况吻合。

针对问题二：考虑风向风速对污染气体扩散过程的影响时，在基于对问题一求解的基础上，在模型一的扩散方程模型中加入风向风速的平流项，由此得到有风情况下的模型，其中分别为风速在x, y方向的分量；对此模型同样采用ADI法求出迭代格式，利用MATLAB编程，求出模型二的数值解，并得到任意时刻污染气体的浓度分布情况。

通过SPSS软件，对附件二所给的原始实际数据与模型二求解得到的模拟值进行显着性检验，检验结果显示该模型与实际情况吻合。

针对问题三：考虑有风时增加高度的影响，此问题即为三维空间的污染气体扩散问题，考虑到三维模型的编程复杂度，而且污染气体的扩散在xoy平面上各向同性，可以将污染气体在y方向的扩散等价为在x方向上的扩散，此时便只需要建立xoz平面上的扩散模型。

在基于对问题二求解的基础上，在模型二的扩散方程中增加高度项，由此得到模型三为，其中为z 方向的扩散系数；对该扩散方程同样采用ADI法求出迭代格式，利用MATLAB编程，求出模型二的数值解，并得到任意时刻污染气体的浓度分布情况。

关键词：污染气体扩散方程ADI法数值解一、问题重述目前，治理雾霾是人们最为关心的热点问题之一。

中国社科院发布的《气候变化绿皮书》中提及，雾霾形成的原因里，重工业、车辆尾气、土方施工都榜上有名，其中钢铁生产也是造成雾霾的重要原因之一。

某钢铁生产集团烟囱污染气体的排放对周边地区大气污染的影响非常大，为了提出更好的治理管理措施，需要对其污染气体扩散的特征进行分析。

现在，我们需要在三种情况下考虑污染气体的扩散过程：1.在不考虑风向和高度影响的情况下，建立模型，模拟某钢铁生产集团的烟囱排放污染气体的扩散过程，假设烟囱的排放速度为300m/s。

2.考虑风向为东北风，平均风速0.6m/s的情况下，模拟污染气体的传播扩散过程。

3.在考虑风向的基础上增加高度的影响，建立模型，模拟污染气体的传播扩散过程。

4.基于上述模型结论，给该钢铁生产集团提供一个污染气体治理建议报告。

二、问题分析钢铁生产集团烟囱污染气体的排放对周边地区大气污染的影响非常大，为了提出更好的治理管理措施，需要对其污染气体扩散的特征进行分析。

污染气体扩散是与时间、空间相关的连续性问题，本文利用微分方程建立模型，通过ADI法（Alternating direction implicit，交替方向隐式法）求解微分方程的迭代式，利用MATLAB编程求其数值解来模拟扩散过程。

对于问题一，在不考虑风速和高度影响的情况下，污染气体的扩散问题即为二维平面点源的扩散问题；考虑到烟囱的排放速度为300m/s，需要在二维扩散方程中加入有源项，由此建立扩散模型。

为检验模型的准确性，需利用SPSS软件，将扩散方程的数值解与实际数据进行显着性分析。

对于问题二，在问题一的基础上，考虑了风向风速的影响，需要在模型一中加入风向风速对应的平流项，由此建立模型二。

为检验模型的准确性，同样需要利用SPSS软件，将模型二的数值解与实际数据进行显着性分析。

对于问题三，基于对问题一、二的建模思路，首先建立不考虑风速风向只考虑高度时的模型，此时应分析扩散系数与高度的关系，建立模型三；然后在模型三的基础上考虑风速风向的影响，此时应分析风速与高度的关系，建立模型四。

为使模型更能符合实际情况下污染气体的扩散过程，本问题还考虑了大气静力稳定度对扩散过程的影响。

综上所述，本问题可以看成是求解偏微分方程中扩散方程的数值解问题。

三、模型假设与符号说明1.假设烟囱的排放速度恒为300m/s；2.假设污染气体中只含有同一类污染物，污染气体的扩散系数相同；3.假设不考虑气体受到的重力和浮力；4.假设在整个扩散过程中污染气体不发生沉降、分解，不发生化学反应；5.假设地面以及地标地物对气体无吸收；四、对问题一的分析与建模4.1问题分析问题一在不考虑风速和高度影响下，烟囱以恒定速度排放污染气体，在无限大平面上，将烟囱口看成是点源，此问题即为二维平面点源的扩散问题，建立二维扩散方程模型，通过ADI法求其数值解，即可得到任意时刻污染气体的浓度分布。

4.2模型一的建立本文将烟囱口当作二维平面上的点源来考虑，分析附件allresults_1所给的数据，可以得出的结论有：污染气体扩散的平面图为1000×1000（单位：m）的正方形区域；以点（500,500）为中心，半径为1m的领域内污染气体的浓度最高，据此确定烟囱口的位置坐标为（500,500）。

假设污染源为连续点源，建立有源项的二维扩散方程：（1）其中，为扩散系数，为时间，为污染源排放污染气体的速度，为浓度对时间的偏导，为浓度对x的二阶偏导，为浓度对y的二阶偏导。

现在考虑初边值条件，初始时刻，平面上各点处浓度均为0，即有边界条件取延拓边界条件，即有根据题意知，在污染源处，污染物排放速度为300m/s，假设烟囱口是以点（500,500）为中心，1m 为半径的圆形区域，在1000×1000的平面内，上述圆形区域任然可以当作点源来处理，在圆形区域内有源项为300，在圆形区域外有源项为0，即有由此得到在不考虑风向和高度影响时连续源的扩散模型（2）另外，根据《中国大百科全书》大气扩散卷，选定扩散系数。

4.3模型一求解本文采用ADI法求解（1）式数值解，将（1）式转化为ADI格式为（a）（b）对（a）式进行转换，令：并将带入（a）式中，得：(3)将(3)式写成矩阵方程组形式：同样对（b）式进行转换，令，并将带入（a）式中，得：(4)将(4)式写成矩阵方程组形式：根据（3）（4）的迭代格式及（2）中初边值条件，运用MATLAB 编程，可以得到任意时刻的污染气体浓度分布，其中截取60s、600s、1800s和3600s时浓度分布如图1.图1（a）60s 时浓度分布图图1（b）600s时浓度分布图图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）两幅图可以看出：由模型一得到的浓度分布图与用实际数据画出的浓度分布图在浓度大小和浓度分布形状上近似一致。

为了更进一步的分析两者之间的误差，现在运用SPSS进一步做模型检验。

取120s 和3600s 时由模型一得到的浓度分布数据与实际数据进行显着性检验。

首先画出模拟值与实际数据的比较图如下：图3（a）120s 模拟值与实际数据比较图图3（b）3600s 时模拟值与实际数据比较图图 3 不同时刻模拟值与实际数据比较图从图3可以直观的看出，模拟值与实际数据之间成明显的线性关系，现在进一步对模拟值与实际数据做线性回归分析，取120s时的模拟数据和实际数据做线性回归分析，SPSS运行结果如下：拟值与实际数据线性回归显着。

且实际值与模拟值之比为1.260，由此可得出模拟值与实际数据基本吻合，即可以用模型一来模拟无风条件下污染气体的扩散过程。

五、对问题二的分析与建模5.1问题分析问题二在问题一的基础上增加了风速和风向的影响，需要在模型一中加入风向风速对应的平流项，同样通过ADI法求其数值解，得到任意时刻污染气体的浓度分布。

5.2模型二的建立假设只考虑水平风向，且风向恒为东北风，风速恒为0.6m/s，将东北风向分解到x ,y 方向上，有其中，分别为x ,y 方向上风速分量。

在模型一的基础上加入风速风向对应的平流项，即有（5）其中，为扩散系数，为时间，为污染源排放污染气体的速度，为浓度对时间的偏导，为浓度对x的偏导，为浓度对y的偏导，为浓度对x的二阶偏导，为浓度对y的二阶偏导。

边界条件取延拓边界条件，即有有源项与问题一中相同，即由此得到在考虑风向影响时连续源的扩散模型（6）5.3模型二的求解运用ADI法求模型二的数值解，ADI法的格式为：令，并将带入（a）式中,对（a）式进行变形，得：（7）将（7）式转化为矩阵方程组形式k=1,2,3,...,K令，并将带入（a）式中,对（b）式进行变形，得：（8）同样，可以将（8）式转化为矩阵方程组形式J=1,2,3,...J根据（7）（8）的迭代格式及（2）中初边值条件，运用MATLAB 编程，可以得到任意时刻的污染气体浓度分布，其中截取60s、600s、1800s和3600s时浓度分布如图4.图4（a）60s 时浓度分布图图4（b）600s 时浓度分布图图4（c）1800s 时浓度分布图图4（d）3600s 时浓度分布图图 4 不同时刻浓度分布图从图4不同时刻浓度分布图可以看出：开始时（t=60s），烟囱口处污染气体的浓度近似为700个单位，远远高于周围其他区域，整个浓度分布区域近似为一条垂直于xoy 平面的直线；随着时间的推移，到t=600s时，烟囱口处污染气体的浓度近似为900个单位，且在（即下风口）区域内的浓度有所增大，整个浓度分布区域由60s的直线状变为半锥形；当t =1800s和t =3600s时，烟囱口处污染气体的浓度任然近似为900个单位，但半锥形底部半圆的面积在不断扩大，到3600s 时已扩散到整个的区域内。