空气污染数学建模.docx

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 时已扩散到整个的区域内。

由此得出的结论是：由于风向风速的存在，烟囱排放的污染气体在平面中不是均匀分布的，在烟囱口处浓度最高，且烟囱口处污染气体浓度随时间增加不明显，在下风口污染气体浓度随时间增加而增大，且污染范围随时间增加而扩大。

5.4模型二的检验

首先通过MATLAB画图，比较模型一的浓度分布图与用实际数据画出的浓度分布图的区别，截取120s 和3600s 的浓度对比图如下：

图 5（a）120s 时模拟值与实际数据浓度分布对比图图 5（b）3600s 时模拟值与实际数据浓度分布对比图

图 5 不同时刻模拟值与实际数据浓度分布对比图

由图5（a）（b）两幅图可以看出：由模型二得到的浓度分布图与用实际数据画出的浓度分布图在浓度大小和浓度分布形状上近似一致。为了更进一步的分析两者之间的误差，现在运用SPSS进一步做模型检验。取120s 和3600s 时由模型二得到的浓度分布数据与实际数据进行显着性检验。首先画出模拟值与实际数据的比较图如下：

图6（a）120s 模拟值与实际数据比较图图6（b）3600s 模拟值与实际数据比较图

图 6 不同时刻模拟值与实际数据比较图

从图6可以直观的看出，模拟值与实际数据之间成明显的线性关系，现在进一步对模拟值与实际数据做线性回归分析，取120s时的模拟数据和实际数据做线性回归分析，SPSS运行结果如下：

拟值与实际数据线性回归显着。且实际值与模拟值之比为1.281，由此可得出模拟值与实际数据基本吻合，即可以用模型二来模拟有风条件下污染气体的扩散过程。

六、对问题三的分析与建模

6.1问题分析

问题三在问题二的基础上增加了高度的影响，假设烟囱高度为10米，此时二维平面模型已不再适用，需要建立三维立体模型来模拟污染气体的扩散过程。考虑到三维模型的编程复杂度，而且污染气体的扩散在xoy平面上各向同性，可以将污染气体在y方向的扩散等价为在x方向上的扩散，此时便只需要建立xoz平面模型。如图7所示。

图7 污染气体扩散的等效二维示意图

对该模型的求解可以采用与问题二同样的分析方法，但要考虑到风速是随高度变化的，需要建立风速与高度的关系式，并将其带入到模型中。考虑到风速为0.6m/s ，[1]中给出的对应的大气稳定度为E,对应的风速随高度变化的关系式为：

（9）其中，h 为高度，为地面风速。在本题中，将地面风速看成0.6m/s .

6.2模型三的建立

基于对问题二的建模过程以及对问题三的分析过程，现建立模型三如下：

（10）其中，为x, y方向的扩散系数，为z方向的扩散系数，为时间，为污染源排放污染气体的速度，为浓度对时间的偏导，为浓度对x方向的偏导，为浓度对y方向的偏导，为对z方向的偏导，为浓度对x的二阶偏导，为浓度对y的二阶偏导。

分别表示风速在x, y方向的分量。具体表达式为

（11）为了简化问题，在求解过程中将z 方向的扩散系数同样当作常数处理，即

由于污染气体在y方向的扩散等价为在x方向上的扩散，（10）式可以写成

（12）边界条件同样取延拓边界条件，即有

假设烟囱高度为10米，于是有源项为

由此得到在考虑风向和高度影响时连续源的扩散模型

（13）6.3模型三的求解

运用ADI法求模型三的数值解，ADI法的格式为：

根据（a）（b）的迭代格式及（13）中初边值条件，运用MATLAB 编程，可以得到任意时刻的污染气体浓度分布，其中截取60s、600s、1800s和3600s时浓度分布如图8.

图 8（a）60s 时浓度分布图图 8（b）600s 时浓度分布图

图 8（c）1800s 时浓度分布图图 8（d）3600s 时浓度分布图

图8 不同时刻浓度分布图

从图8不同时刻浓度分布图可以看出：任意时刻，烟囱口附近的污染气体浓度最高；在x 方向，以烟囱口为中心，污染气体的浓度在左右两边呈不对称分布。开始时（t=60s），污染气体在左边扩散到200米处，在右边扩散到580米处；随着时间的推移，到t=600s 时，污染气体在左边扩散到100米处，在右边仍为580米处；当t =1800s时，污染气体在左边扩散到50米处，在右边没有什么变化，仍为580米处；当t =3600s时，污染气体在左边扩散至区域的边缘处，右边仍然在580米处。在z 方向，以烟囱口为中心，污染气体的浓度在上下两侧同样呈不对称分布，任意时刻，下侧的浓度扩散范围始终大于上侧，到3600s 时，污染气体已扩散到距地面2m 左右。

由此得出的结论是：由于风向风速及烟囱高度的影响，烟囱排放的污染气体在空间中不是均匀分布的，在烟囱口处浓度最高，随着时间的推移，上风口处污染气体的浓度及污染范围变化不大，而在下风口处，污染气体的浓度和污染范围在不断地增大，到3600s 时，污染气体已扩散至近地面处。

现在将烟囱的高度设为20米，在其他条件不变的情况下分析污染气体的扩散过程，并将其与烟囱高度为10米情况下进行比较，截取3600s 时比较结果如下图：

图 9（a）烟囱高度为10米的扩散结果图 9（b）烟囱高度为20米的扩散结果

图9 不同烟囱高度的扩散结果

从图9的a ,b两幅图可以看出：当烟囱高度为20米时，污染气体在z 方向能向下扩散5米左右，即地表及地表以上15米的范围内未受到污染气体的影响。

由此得出结论：增加烟囱高度能使污染气体不扩散至近地面，而是在十几米的低空中因大气运动扩散出去，不会对地面人群造成影响。

七、给钢铁生产集团治理污染气体的建议报告

目前，治理雾霾是人们最为关心的热点问题之一，雾霾形成的重要原因之一便是钢铁生产，为了减小钢铁生产集团烟囱污染气体的排放对周边地区大气污染的影响，基于污染气体扩散模型，现提出以下建议：

1.减小污染气体排放速度。当排放量一定时，减小污染气体排放速度，在单位时间内大

气中污染物浓度的增量较小，减少了污染气体在环境大气中的积累，利于污染物的扩散，使得大气在较短时间内恢复自净能力。

2.考虑选址问题。钢铁厂适宜建在城市的下风口，因此在建厂之前需要调查了解城市的

盛行风向，将厂址选在远离居民区的下风向，以此减小污染气体对居民及周围大气的影响。

3.合理增加烟囱高度。在满足国家有关烟囱高度标准的前提下，适当增加烟囱高度，可

以使得污染气体的污染范围保持在低空，不至于扩散到近地面，从而保证了地面人群不受污染气体的影响。

污染气体的治理是一个长期的过程，钢铁生产集团在基于以上建议的基础上，制订出更好的治理管理措施，最终使得钢铁生产和环境治理达到双赢的结果。

八、模型优缺点分析

8.1模型优点

（1）问题一中建立了无风及不考虑高度时的污染气体扩散模型。模型一建立过程中，考虑了污染气体排放速度的实际情况，在二维扩散方程中加入了有源项，对

该二维扩散方程的求解采用ADI法，得到迭代关系式，通过MATLAB编程，求

出二维扩散方程的数值解，并得到任意时刻污染气体的浓度分布情况。通过附件

一数据的检验，该模型能够很好的描述污染气体的扩散过程。

（2）问题二在问题一的基础上考虑了风向风速的影响，在模型一的二维扩散方程中加入了风向风速的平流项，得到模型二。同样采用ADI法，得到迭代关系式，利用MATLAB编程，求出模型二的数值解，并得到任意时刻污染气体的浓度分

布情况。通过附件二数据的检验，该模型能够很好的描述有风情况下污染气体的

扩散过程。

（3）问题三在问题二的基础上考虑了高度的影响，在模型三的建立过程中，将复杂的三维空间扩散模型转化为二维扩散模型，简化了编程的复杂度，根据模型三，利用MATLAB编程，模拟出了高度对污染气体扩散过程的影响。

8.2模型缺点

（1）在本文中将污染气体扩散系数当作常数处理，未考虑扩散系数与高度等因素的关系，导致模型与实际情况出现偏差。

（2）在本文中只考虑了同一类污染气体的扩散过程，未考虑多种污染气体同时存在时的扩散情况，导致模型与实际出现偏差。

九、模型推广与应用

模型一建立了无风及不考虑高度时的污染气体扩散模型，该模型考虑了污染气体排放速度的实际情况，在二维扩散方程中加入了有源项，采用ADI法得到扩散方程的迭代式，利用MATLAB编程，求出其数值解，得到任意时刻污染气体的浓度分布。经检验该模型与实际相符，所以该模型可用于求解平面上无风情况下的污染气体扩散过程。

模型二建立了有风情况下的污染气体扩散模型，采用与模型一同样的处理方法，得到任意时刻污染气体的浓度分布。经检验该模型与实际相符，所以该模型可用于求解平面上有风情况下的污染气体扩散过程。

模型三综合考虑了风向风速及高度等因素的影响，使模型能够更好的运用于实际情况中，能较好的模拟真实情况下污染气体的扩散过程，为钢铁生产集团提供理论依据。

参考文献

[1]王志春宋丽莉何秋生刘爱君刘荣叶燕翔，风速随高度变化的曲线模型分析，热

带气象学报，第23卷第6期：690-692页，2007

[2]姜启源谢金星叶俊，数学模型（第四版），北京：高等教育出版社，2011

[3]Rafael B. Storch Luiz C.G. Pimentel Helcio R.B. Orlande，Identification of

atmospheric boundary layer parameters by inverse problem , atmospheric environment 41 :1417-1425,2007

附录：

1.无风模型

clc;

clear;

% N is the total time steps

% J is the total mesh points in X

% K is the total mesh points in Y

% for space intervals

N=60;

xb=0.0;

xe=1000.0;

yb=0.0;

ye=1000.0;

J=101;

K=101;

% for time intervals

t0=0;

dt=1;

tf=dt*N;

% the number of mesh points

J01=J-1;

K01=K-1;

dx=(xe-xb)/(J-1);

dy=(ye-yb)/(K-1);

%dt=(tf-t0)/(N-1);

a=10;

nux=dt*a/(dx^2);

nuy=dt*a/(dy^2);

%to obtain the coordinates

x=xb:dx:xe;

y=yb:dy:ye;

u1=zeros(J,K);

u2=zeros(J,K);

%Initial condition:

t=0;

for k=1:K

for j=1:J

u1(j,k)=0;

end

end

%-------Big loop for time t -------

for n =1:N

t1 = t + dt;t2 = t + dt/2;

%-----sweep in x-direction ---------

%boundary conditions

u2(:,1)=0;

u2(:,K)=0;

for k = 2:K01, %look for fixed y(k)

A = sparse(J,J);b = zeros(J,1);

for j = 2:J01,

b(j) = (u1(j,k-1)-2*u1(j,k)+u1(j,k+1))*nuy/2+f(x(j),y(k))*dt/2+u1(j,k);

A(j,j-1) = -nux/2;

A(j,j) = 1+nux;

A(j,j+1) = -nux/2;

end

A(1,1)=1;

A(J,J)=1;

b(1)=0;%these are for zero boundary conditions

b(J)=0;

ut = A\b;%solve the tri-diagonal matrix.

u2(:,k)=ut;

end %finishi x-sweep

%-----sweep in y-direction ---------

%boundary conditions

u1(1,:)=0;

u1(J,:)=0;

for j = 2:J01, %look for fixed x(k)

A = sparse(K,K);b = zeros(K,1);

for k = 2:K01,

b(k) = (u2(j-1,k)-2*u2(j,k)+u2(j+1,k))*nux/2+f(x(j),y(k))*dt/2+u2(j,k);

A(k,k-1) = -nuy/2;

A(k,k) = 1+nuy;

A(k,k+1) = -nuy/2;

end

A(1,1)=1;

A(K,K)=1;

b(1)=0;%these are for zero boundary conditions

b(K)=0;

ut = A\b;%solve the tri-diagonal matrix.

u1(j,:)=ut;

end %finishi x-sweep

end

mesh(x,y,u1)

subplot(1,2,1)

mesh(x,y,u1)

xlabel('x')

ylabel('y')

zlabel('浓度c')

title('模型求解图')

load ('a160s.txt');

x1=a160s(:,1);

y1=a160s(:,2);

c1=a160s(:,3);

subplot(1,2,2)

plot3(x1,y1,c1,'b')

xlabel('x')

ylabel('y')

zlabel('浓度c')

title('原始数据曲线图')

suptitle('t=60s')

2.无风时源项f

function [ result ] = f(x,y)

if (x-500)^2+(y-500)^2<=1

result=300;

else

result=0;

end

end

3.有风模型

clc;

clear;

% N is the total time steps

% J is the total mesh points in X

% K is the total mesh points in Y

% for space intervals

N=360;

xb=0.0;

xe=1000.0;

yb=0.0;

ye=1000.0;

J=101;

K=101;

% for time intervals

t0=0;

dt=1;

tf=dt*N;

% the number of mesh points

J01=J-1;

K01=K-1;

dx=(xe-xb)/(J-1);

dy=(ye-yb)/(K-1);

%dt=(tf-t0)/(N-1);

a=10;

m1=sqrt(2)*0.3;

m2=sqrt(2)*0.3;

nux=dt*a/(dx^2);

nuy=dt*a/(dy^2);

mux=dt*m1/(2*dx);

muy=dt*m2/(2*dy);

%to obtain the coordinates

x=xb:dx:xe;

y=yb:dy:ye;

u1=zeros(J,K);

u2=zeros(J,K);

%Initial condition:

t=0;

for k=1:K

for j=1:J

u1(j,k)=0;

end

end

%-------Big loop for time t -------

for n =1:N

t1 = t + dt;t2 = t + dt/2;

%-----sweep in x-direction ---------

%boundary conditions

u2(:,1)=0;

u2(:,K)=0;

for k = 2:K01, %look for fixed y(k)

A = sparse(J,J);b = zeros(J,1);

for j = 2:J01,

b(j) = (u1(j,k-1)-2*u1(j,k)+u1(j,k+1))*nuy/2+f(x(j),y(k))*dt/2+u1(j,k)+muy/2*(u1(j,k+1)-u1(j,k-1) );

A(j,j-1) = -nux/2+mux/2;

A(j,j) = 1+nux;

A(j,j+1) = -nux/2-mux/2;

end

A(1,1)=1;

A(J,J)=1;

b(1)=0;%these are for zero boundary conditions

b(J)=0;

ut = A\b;%solve the tri-diagonal matrix.

u2(:,k)=ut;

end %finishi x-sweep

%-----sweep in y-direction ---------

%boundary conditions

for j = 2:J01, %look for fixed x(k)

A = sparse(K,K);b = zeros(K,1);

for k = 2:K01,

b(k) = (u2(j-1,k)-2*u2(j,k)+u2(j+1,k))*nux/2+f(x(j),y(k))*dt/2+u2(j,k)+mux/2*(u2(j+1,k)-u2(j-1,k) );

A(k,k-1) = -nuy/2+muy/2;

A(k,k) = 1+nuy;

A(k,k+1) = -nuy/2-muy/2;

end

A(1,1)=1;

A(K,K)=1;

b(1)=0;%these are for zero boundary conditions

b(K)=0;

ut = A\b;%solve the tri-diagonal matrix.

u1(j,:)=ut;

end %finishi x-sweep

end

mesh(x,y,u1)

4.有风时的源项f2

function [ result ] = f2(x,z)

if (x-500)^2+(z-10)^2<=1

result=300;

else

result=0;

end

End

5.原始数据浓度的导出

clc;

clear;

load ('a180.txt');

c1=a180(:,3);

dlmwrite('a1801s.txt',c1)

6.考虑风速风向和高度的模型

clc;

clear;

% N is the total time steps

% J is the total mesh points in X

% K is the total mesh points in Y

% for space intervals

N=60;

xb=0.0;

xe=1000.0;

J=101;

K=101;

% for time intervals

t0=0;

dt=1;

tf=dt*N;

% the number of mesh points

J01=J-1;

K01=K-1;

dx=(xe-xb)/(J-1);

dz=(ze-zb)/(K-1);

%dt=(tf-t0)/(N-1);

a1=20;

a2=10;

nux=dt*a1/(dx^2);

nuz=dt*a2/(dz^2);

%to obtain the coordinates

x=xb:dx:xe;

z=zb:dz:ze;

u1=zeros(J,K);

u2=zeros(J,K);

%Initial condition:

t=0;

for k=1:K

for j=1:J

u1(j,k)=0;

end

end

%-------Big loop for time t -------

for n =1:N

t1 = t + dt;t2 = t + dt/2;

%-----sweep in x-direction ---------

%boundary conditions

u2(:,1)=0;

u2(:,K)=0;

for k = 2:K01, %look for fixed y(k)

m=sqrt(2)*0.3*(k*dz/10)^0.25;

mux=dt*m/(2*dx);

muz=dt*m/(2*dz);

A = sparse(J,J);b = zeros(J,1);

for j = 2:J01,

b(j) = (u1(j,k-1)-2*u1(j,k)+u1(j,k+1))*nuz/2+f2(x(j),z(k))*dt/2+u1(j,k)-muz/2*(u1(j,k+1)-u1(j,k-1)

);

A(j,j-1) = -nux/2-mux/2;

A(j,j) = 1+nux;

A(j,j+1) = -nux/2+mux/2;

end

A(1,1)=1;

A(J,J)=1;

b(1)=0;%these are for zero boundary conditions

b(J)=0;

ut = A\b;%solve the tri-diagonal matrix.

u2(:,k)=ut;

end %finishi x-sweep

%-----sweep in y-direction ---------

%boundary conditions

u1(1,:)=0;

u1(J,:)=0;

for j = 2:J01, %look for fixed x(k)

A = sparse(K,K);b = zeros(K,1);

for k = 2:K01

m=sqrt(2)*0.3*(k*dz/10)^0.25;

mux=dt*m/(2*dx);

muz=dt*m/(2*dz);

b(k) = (u2(j-1,k)-2*u2(j,k)+u2(j+1,k))*nux/2+f2(x(j),z(k))*dt/2+u2(j,k)-mux/2*(u2(j+1,k)-u2(j-1,k) );

A(k,k-1) = -nuz/2-muz/2;

A(k,k) = 1+nuz;

A(k,k+1) = -nuz/2+muz/2;

end

A(1,1)=1;

A(K,K)=1;

b(1)=0;%these are for zero boundary conditions

b(K)=0;

ut = A\b;%solve the tri-diagonal matrix.

u1(j,:)=ut;

end %finishi x-sweep

end

contour(x,z,u1)

xlabel('x')

ylabel('z')

zlabel('浓度c')

title('模型求解图')