空氣污染數(shù)學(xué)建模27頁(yè)

1、A.污染氣體的傳播擴(kuò)散摘要鋼鐵生產(chǎn)排放的污染氣體是造成霧霾的重要原因之一，研究污染氣體的擴(kuò)散特征，正確模擬污染氣體的擴(kuò)散過(guò)程，能夠?yàn)殇撹F生產(chǎn)集團(tuán)提出更好的治理管理措施，具有實(shí)際意義。針對(duì)問(wèn)題一：污染氣體的排放速度為300m/s,在不考慮風(fēng)向風(fēng)速及高度影響的情況下，此問(wèn)題即為二維平面的連續(xù)點(diǎn)源擴(kuò)散問(wèn)題，由此在二維xoy 平面上建立連續(xù)點(diǎn)源擴(kuò)散方程模型，其中 為氣體擴(kuò)散系數(shù)，本文中取為常數(shù)10，f(x,y,t ) 為污染氣體的排放速度，在本文中恒為300m/s ；對(duì)上述偏微分方程模型，本文采用ADI法（Alternating direction implicit，交替方向隱式法）求解出迭代格式，利

2、用MATLAB編程，求出模型一的數(shù)值解，并得到任意時(shí)刻污染氣體的濃度分布情況。通過(guò)SPSS軟件，對(duì)附件一所給的原始實(shí)際數(shù)據(jù)與模型一求解得到的模擬值進(jìn)行顯著性檢驗(yàn)，檢驗(yàn)結(jié)果顯示該模型與實(shí)際情況吻合。針對(duì)問(wèn)題二：考慮風(fēng)向風(fēng)速對(duì)污染氣體擴(kuò)散過(guò)程的影響時(shí)，在基于對(duì)問(wèn)題一求解的基礎(chǔ)上，在模型一的擴(kuò)散方程模型中加入風(fēng)向風(fēng)速的平流項(xiàng)，由此得到有風(fēng)情況下的模型，其中 分別為風(fēng)速在x, y方向的分量；對(duì)此模型同樣采用ADI法求出迭代格式，利用MATLAB編程，求出模型二的數(shù)值解，并得到任意時(shí)刻污染氣體的濃度分布情況。通過(guò)SPSS軟件，對(duì)附件二所給的原始實(shí)際數(shù)據(jù)與模型二求解得到的模擬值進(jìn)行顯著性檢驗(yàn)，檢驗(yàn)結(jié)果顯示

3、該模型與實(shí)際情況吻合。針對(duì)問(wèn)題三：考慮有風(fēng)時(shí)增加高度的影響，此問(wèn)題即為三維空間的污染氣體擴(kuò)散問(wèn)題，考慮到三維模型的編程復(fù)雜度，而且污染氣體的擴(kuò)散在xoy平面上各向同性，可以將污染氣體在y方向的擴(kuò)散等價(jià)為在x方向上的擴(kuò)散，此時(shí)便只需要建立xoz平面上的擴(kuò)散模型。在基于對(duì)問(wèn)題二求解的基礎(chǔ)上，在模型二的擴(kuò)散方程中增加高度項(xiàng)，由此得到模型三為，其中 為z 方向的擴(kuò)散系數(shù)；對(duì)該擴(kuò)散方程同樣采用ADI法求出迭代格式，利用MATLAB編程，求出模型二的數(shù)值解，并得到任意時(shí)刻污染氣體的濃度分布情況。關(guān)鍵詞：污染氣體 擴(kuò)散方程 ADI法 數(shù)值解1、 問(wèn)題重述目前，治理霧霾是人們最為關(guān)心的熱點(diǎn)問(wèn)題之一。中國(guó)社科院

4、發(fā)布的氣候變化綠皮書中提及，霧霾形成的原因里，重工業(yè)、車輛尾氣、土方施工都榜上有名，其中鋼鐵生產(chǎn)也是造成霧霾的重要原因之一。某鋼鐵生產(chǎn)集團(tuán)煙囪污染氣體的排放對(duì)周邊地區(qū)大氣污染的影響非常大，為了提出更好的治理管理措施，需要對(duì)其污染氣體擴(kuò)散的特征進(jìn)行分析?，F(xiàn)在，我們需要在三種情況下考慮污染氣體的擴(kuò)散過(guò)程：1. 在不考慮風(fēng)向和高度影響的情況下，建立模型，模擬某鋼鐵生產(chǎn)集團(tuán)的煙囪排放污染氣體的擴(kuò)散過(guò)程，假設(shè)煙囪的排放速度為300m/s。2. 考慮風(fēng)向?yàn)闁|北風(fēng)，平均風(fēng)速0.6m/s的情況下，模擬污染氣體的傳播擴(kuò)散過(guò)程。3. 在考慮風(fēng)向的基礎(chǔ)上增加高度的影響，建立模型，模擬污染氣體的傳播擴(kuò)散過(guò)程。4. 基

5、于上述模型結(jié)論，給該鋼鐵生產(chǎn)集團(tuán)提供一個(gè)污染氣體治理建議報(bào)告。2、 問(wèn)題分析鋼鐵生產(chǎn)集團(tuán)煙囪污染氣體的排放對(duì)周邊地區(qū)大氣污染的影響非常大，為了提出更好的治理管理措施，需要對(duì)其污染氣體擴(kuò)散的特征進(jìn)行分析。污染氣體擴(kuò)散是與時(shí)間、空間相關(guān)的連續(xù)性問(wèn)題，本文利用微分方程建立模型，通過(guò)ADI法（Alternating direction implicit，交替方向隱式法）求解微分方程的迭代式，利用MATLAB編程求其數(shù)值解來(lái)模擬擴(kuò)散過(guò)程。對(duì)于問(wèn)題一，在不考慮風(fēng)速和高度影響的情況下，污染氣體的擴(kuò)散問(wèn)題即為二維平面點(diǎn)源的擴(kuò)散問(wèn)題；考慮到煙囪的排放速度為300m/s，需要在二維擴(kuò)散方程中加入有源項(xiàng)，由此建立擴(kuò)

6、散模型。為檢驗(yàn)?zāi)Ｐ偷臏?zhǔn)確性，需利用SPSS軟件，將擴(kuò)散方程的數(shù)值解與實(shí)際數(shù)據(jù)進(jìn)行顯著性分析。對(duì)于問(wèn)題二，在問(wèn)題一的基礎(chǔ)上，考慮了風(fēng)向風(fēng)速的影響，需要在模型一中加入風(fēng)向風(fēng)速對(duì)應(yīng)的平流項(xiàng)，由此建立模型二。為檢驗(yàn)?zāi)Ｐ偷臏?zhǔn)確性，同樣需要利用SPSS軟件，將模型二的數(shù)值解與實(shí)際數(shù)據(jù)進(jìn)行顯著性分析。對(duì)于問(wèn)題三，基于對(duì)問(wèn)題一、二的建模思路，首先建立不考慮風(fēng)速風(fēng)向只考慮高度時(shí)的模型，此時(shí)應(yīng)分析擴(kuò)散系數(shù)與高度的關(guān)系，建立模型三；然后在模型三的基礎(chǔ)上考慮風(fēng)速風(fēng)向的影響，此時(shí)應(yīng)分析風(fēng)速與高度的關(guān)系，建立模型四。為使模型更能符合實(shí)際情況下污染氣體的擴(kuò)散過(guò)程，本問(wèn)題還考慮了大氣靜力穩(wěn)定度對(duì)擴(kuò)散過(guò)程的影響。綜上所述，本問(wèn)

7、題可以看成是求解偏微分方程中擴(kuò)散方程的數(shù)值解問(wèn)題。3、 模型假設(shè)與符號(hào)說(shuō)明1. 假設(shè)煙囪的排放速度恒為300m/s；2. 假設(shè)污染氣體中只含有同一類污染物，污染氣體的擴(kuò)散系數(shù)相同；3. 假設(shè)不考慮氣體受到的重力和浮力；4. 假設(shè)在整個(gè)擴(kuò)散過(guò)程中污染氣體不發(fā)生沉降、分解，不發(fā)生化學(xué)反應(yīng)；5. 假設(shè)地面以及地標(biāo)地物對(duì)氣體無(wú)吸收；6. 假設(shè)在有風(fēng)情況下，風(fēng)向水平，風(fēng)速風(fēng)向恒定；u : 污染物濃度f(wàn) : 污染物排放速度v : 風(fēng)速 : 水平方向擴(kuò)散系數(shù)： 風(fēng)速在x 方向的分量： 風(fēng)速在y 方向的分量 ： 垂直方向擴(kuò)散系數(shù)4、 對(duì)問(wèn)題一的分析與建模4.1 問(wèn)題分析問(wèn)題一在不考慮風(fēng)速和高度影響下，煙囪以恒

8、定速度排放污染氣體，在無(wú)限大平面上，將煙囪口看成是點(diǎn)源，此問(wèn)題即為二維平面點(diǎn)源的擴(kuò)散問(wèn)題，建立二維擴(kuò)散方程模型，通過(guò)ADI法求其數(shù)值解，即可得到任意時(shí)刻污染氣體的濃度分布。4.2 模型一的建立本文將煙囪口當(dāng)作二維平面上的點(diǎn)源來(lái)考慮，分析附件allresults_1所給的數(shù)據(jù)，可以得出的結(jié)論有：污染氣體擴(kuò)散的平面圖為10001000（單位：m）的正方形區(qū)域；以點(diǎn)（500,500）為中心，半徑為1m的領(lǐng)域內(nèi)污染氣體的濃度最高，據(jù)此確定煙囪口的位置坐標(biāo)為（500,500）。假設(shè)污染源為連續(xù)點(diǎn)源，建立有源項(xiàng)的二維擴(kuò)散方程： （1）其中，為擴(kuò)散系數(shù)， 為時(shí)間，為污染源排放污染氣體的速度，為濃度對(duì)時(shí)間的偏

9、導(dǎo)， 為濃度對(duì)x的二階偏導(dǎo)，為濃度對(duì)y的二階偏導(dǎo)。現(xiàn)在考慮初邊值條件，初始時(shí)刻，平面上各點(diǎn)處濃度均為0，即有邊界條件取延拓邊界條件，即有 根據(jù)題意知，在污染源處，污染物排放速度為300m/s，假設(shè)煙囪口是以點(diǎn)（500,500）為中心，1m 為半徑的圓形區(qū)域，在10001000的平面內(nèi)，上述圓形區(qū)域任然可以當(dāng)作點(diǎn)源來(lái)處理，在圓形區(qū)域內(nèi)有源項(xiàng)為300，在圓形區(qū)域外有源項(xiàng)為0，即有由此得到在不考慮風(fēng)向和高度影響時(shí)連續(xù)源的擴(kuò)散模型 （2）另外，根據(jù)中國(guó)大百科全書大氣擴(kuò)散卷，選定擴(kuò)散系數(shù)。4.3 模型一求解本文采用ADI法求解（1）式數(shù)值解，將（1）式轉(zhuǎn)化為ADI格式為 （a）（b）對(duì)（a）式進(jìn)行轉(zhuǎn)換，

10、令： 并將帶入（a）式中，得： (3)將(3)式寫成矩陣方程組形式： 同樣對(duì)（b）式進(jìn)行轉(zhuǎn)換，令，并將帶入（a）式中，得： (4)將(4)式寫成矩陣方程組形式： 根據(jù)（3）（4）的迭代格式及（2）中初邊值條件，運(yùn)用MATLAB 編程，可以得到任意時(shí)刻的污染氣體濃度分布，其中截取60s、600s、1800s和3600s時(shí)濃度分布如圖1.圖1（a）60s 時(shí)濃度分布圖圖 1（b）600s時(shí)濃度分布圖圖 1（c）1800s 時(shí)濃度分布圖圖 1（d）3600s 時(shí)濃度分布圖 1 不同時(shí)刻濃度分布圖從圖1不同時(shí)刻濃度分布圖可以看出：開始時(shí)（t=60s），煙囪口處污染氣體的濃度近似為800個(gè)單位，遠(yuǎn)遠(yuǎn)高于

11、周圍其他區(qū)域，整個(gè)濃度分布區(qū)域近似為一條垂直于xoy平面的直線；隨著時(shí)間的推移，到t=600s時(shí)，煙囪口處污染氣體的濃度近似為1100個(gè)單位，且周圍其他區(qū)域內(nèi)的濃度有所增大，整個(gè)濃度分布區(qū)域由60s的直線狀變?yōu)槌噬霞庀聢A的塔狀；當(dāng)t =1800s時(shí)，煙囪口處污染氣體的濃度近似為1300個(gè)單位；當(dāng)t =3600s時(shí)，煙囪口處污染氣體的濃度已達(dá)到1700個(gè)單位左右，從“塔”的形狀來(lái)看，“塔尖”在不斷地增高，“塔底”在不斷地向外擴(kuò)展，到3600s 時(shí)已擴(kuò)散到半徑為800米的圓形區(qū)域內(nèi)。4.4 模型一的檢驗(yàn)首先通過(guò)MATLAB畫圖，比較模型一的濃度分布圖與用實(shí)際數(shù)據(jù)畫出的濃度分布圖的區(qū)別，截取120s

12、 和3600s 的濃度對(duì)比圖如下：圖 2（a）120s 時(shí)模擬值與實(shí)際數(shù)據(jù)濃度分布對(duì)比圖圖 2（b）3600s 時(shí)模擬值與實(shí)際數(shù)據(jù)濃度分布對(duì)比圖圖 2 不同時(shí)刻模擬值與實(shí)際數(shù)據(jù)濃度分布對(duì)比圖由圖2（a）（b）兩幅圖可以看出：由模型一得到的濃度分布圖與用實(shí)際數(shù)據(jù)畫出的濃度分布圖在濃度大小和濃度分布形狀上近似一致。為了更進(jìn)一步的分析兩者之間的誤差，現(xiàn)在運(yùn)用SPSS進(jìn)一步做模型檢驗(yàn)。取120s 和3600s 時(shí)由模型一得到的濃度分布數(shù)據(jù)與實(shí)際數(shù)據(jù)進(jìn)行顯著性檢驗(yàn)。首先畫出模擬值與實(shí)際數(shù)據(jù)的比較圖如下：圖 3（a）120s 模擬值與實(shí)際數(shù)據(jù)比較圖圖 3（b）3600s 時(shí)模擬值與實(shí)際數(shù)據(jù)比較圖圖 3 不

13、同時(shí)刻模擬值與實(shí)際數(shù)據(jù)比較圖從圖3可以直觀的看出，模擬值與實(shí)際數(shù)據(jù)之間成明顯的線性關(guān)系，現(xiàn)在進(jìn)一步對(duì)模擬值與實(shí)際數(shù)據(jù)做線性回歸分析，取120s時(shí)的模擬數(shù)據(jù)和實(shí)際數(shù)據(jù)做線性回歸分析，SPSS運(yùn)行結(jié)果如下：Anovaa模型平方和df均方FSig.1回歸12176482.538112176482.5382966236.658.000b殘差41867.174101994.105總計(jì)12218349.71110200a. 因變量: V1b. 預(yù)測(cè)變量: (常量), V2。系數(shù)a模型非標(biāo)準(zhǔn)化系數(shù)標(biāo)準(zhǔn)系數(shù)tSig.B 的 95.0% 置信區(qū)間B標(biāo)準(zhǔn) 誤差試用版下限上限實(shí)際值(常量)-.089.020-4.4

14、32.000-.129-.050模擬值1.260.001.9981722.277.0001.2591.262a. 因變量: V1從運(yùn)行結(jié)果中可以看出：回歸檢驗(yàn)的p 值為0.000，小于0.05，拒絕原假設(shè)，即模擬值與實(shí)際數(shù)據(jù)線性回歸顯著。且實(shí)際值與模擬值之比為1.260，由此可得出模擬值與實(shí)際數(shù)據(jù)基本吻合，即可以用模型一來(lái)模擬無(wú)風(fēng)條件下污染氣體的擴(kuò)散過(guò)程。5、 對(duì)問(wèn)題二的分析與建模5.1 問(wèn)題分析問(wèn)題二在問(wèn)題一的基礎(chǔ)上增加了風(fēng)速和風(fēng)向的影響，需要在模型一中加入風(fēng)向風(fēng)速對(duì)應(yīng)的平流項(xiàng)，同樣通過(guò)ADI法求其數(shù)值解，得到任意時(shí)刻污染氣體的濃度分布。5.2 模型二的建立假設(shè)只考慮水平風(fēng)向，且風(fēng)向恒為東北

15、風(fēng)，風(fēng)速恒為0.6m/s，將東北風(fēng)向分解到x ,y 方向上，有其中， 分別為x ,y 方向上風(fēng)速分量。在模型一的基礎(chǔ)上加入風(fēng)速風(fēng)向?qū)?yīng)的平流項(xiàng)，即有 （5）其中，為擴(kuò)散系數(shù)， 為時(shí)間，為污染源排放污染氣體的速度，為濃度對(duì)時(shí)間的偏導(dǎo)，為濃度對(duì)x的偏導(dǎo)，為濃度對(duì)y的偏導(dǎo)，為濃度對(duì)x的二階偏導(dǎo)，為濃度對(duì)y的二階偏導(dǎo)。邊界條件取延拓邊界條件，即有 有源項(xiàng) 與問(wèn)題一中相同，即由此得到在考慮風(fēng)向影響時(shí)連續(xù)源的擴(kuò)散模型 （6）5.3 模型二的求解運(yùn)用ADI法求模型二的數(shù)值解，ADI法的格式為：令 ，并將帶入（a）式中,對(duì)（a）式進(jìn)行變形，得：（7）將（7）式轉(zhuǎn)化為矩陣方程組形式 k=1,2,3,.,K令，并

16、將帶入（a）式中,對(duì)（b）式進(jìn)行變形，得：（8）同樣，可以將（8）式轉(zhuǎn)化為矩陣方程組形式 J=1,2,3,.J根據(jù)（7）（8）的迭代格式及（2）中初邊值條件，運(yùn)用MATLAB 編程，可以得到任意時(shí)刻的污染氣體濃度分布，其中截取60s、600s、1800s和3600s時(shí)濃度分布如圖4.圖 4（a）60s 時(shí)濃度分布圖圖 4（b）600s 時(shí)濃度分布圖圖 4（c）1800s 時(shí)濃度分布圖圖 4（d）3600s 時(shí)濃度分布圖圖 4 不同時(shí)刻濃度分布圖從圖4不同時(shí)刻濃度分布圖可以看出：開始時(shí)（t=60s），煙囪口處污染氣體的濃度近似為700個(gè)單位，遠(yuǎn)遠(yuǎn)高于周圍其他區(qū)域，整個(gè)濃度分布區(qū)域近似為一條垂直于

17、xoy平面的直線；隨著時(shí)間的推移，到t=600s時(shí)，煙囪口處污染氣體的濃度近似為900個(gè)單位，且在 （即下風(fēng)口）區(qū)域內(nèi)的濃度有所增大，整個(gè)濃度分布區(qū)域由60s的直線狀變?yōu)榘脲F形；當(dāng)t =1800s和t =3600s時(shí)，煙囪口處污染氣體的濃度任然近似為900個(gè)單位，但半錐形底部半圓的面積在不斷擴(kuò)大，到3600s 時(shí)已擴(kuò)散到整個(gè)的區(qū)域內(nèi)。由此得出的結(jié)論是：由于風(fēng)向風(fēng)速的存在，煙囪排放的污染氣體在平面中不是均勻分布的，在煙囪口處濃度最高，且煙囪口處污染氣體濃度隨時(shí)間增加不明顯，在下風(fēng)口污染氣體濃度隨時(shí)間增加而增大，且污染范圍隨時(shí)間增加而擴(kuò)大。5.4 模型二的檢驗(yàn)首先通過(guò)MATLAB畫圖，比較模型一的

18、濃度分布圖與用實(shí)際數(shù)據(jù)畫出的濃度分布圖的區(qū)別，截取120s 和3600s 的濃度對(duì)比圖如下：圖 5（a）120s 時(shí)模擬值與實(shí)際數(shù)據(jù)濃度分布對(duì)比圖圖 5（b）3600s 時(shí)模擬值與實(shí)際數(shù)據(jù)濃度分布對(duì)比圖圖 5 不同時(shí)刻模擬值與實(shí)際數(shù)據(jù)濃度分布對(duì)比圖由圖5（a）（b）兩幅圖可以看出：由模型二得到的濃度分布圖與用實(shí)際數(shù)據(jù)畫出的濃度分布圖在濃度大小和濃度分布形狀上近似一致。為了更進(jìn)一步的分析兩者之間的誤差，現(xiàn)在運(yùn)用SPSS進(jìn)一步做模型檢驗(yàn)。取120s 和3600s 時(shí)由模型二得到的濃度分布數(shù)據(jù)與實(shí)際數(shù)據(jù)進(jìn)行顯著性檢驗(yàn)。首先畫出模擬值與實(shí)際數(shù)據(jù)的比較圖如下：圖 6（a）120s 模擬值與實(shí)際數(shù)據(jù)比較圖

19、圖 6（b）3600s 模擬值與實(shí)際數(shù)據(jù)比較圖圖 6 不同時(shí)刻模擬值與實(shí)際數(shù)據(jù)比較圖從圖6可以直觀的看出，模擬值與實(shí)際數(shù)據(jù)之間成明顯的線性關(guān)系，現(xiàn)在進(jìn)一步對(duì)模擬值與實(shí)際數(shù)據(jù)做線性回歸分析，取120s時(shí)的模擬數(shù)據(jù)和實(shí)際數(shù)據(jù)做線性回歸分析，SPSS運(yùn)行結(jié)果如下：Anovaa模型平方和df均方FSig.1回歸5789448.83515789448.8351428052.592.000b殘差41347.629101994.054總計(jì)5830796.46410200a. 因變量: V1b. 預(yù)測(cè)變量: (常量), V2。系數(shù)a模型非標(biāo)準(zhǔn)化系數(shù)標(biāo)準(zhǔn)系數(shù)tSig.B標(biāo)準(zhǔn) 誤差試用版實(shí)際值(常量)-.074.

20、020-3.715.000模擬值1.281.001.9961195.012.000a. 因變量: V1從運(yùn)行結(jié)果中可以看出：回歸檢驗(yàn)的p 值為0.000，小于0.05，拒絕原假設(shè)，即模擬值與實(shí)際數(shù)據(jù)線性回歸顯著。且實(shí)際值與模擬值之比為1.281，由此可得出模擬值與實(shí)際數(shù)據(jù)基本吻合，即可以用模型二來(lái)模擬有風(fēng)條件下污染氣體的擴(kuò)散過(guò)程。6、 對(duì)問(wèn)題三的分析與建模6.1 問(wèn)題分析問(wèn)題三在問(wèn)題二的基礎(chǔ)上增加了高度的影響，假設(shè)煙囪高度為10米，此時(shí)二維平面模型已不再適用，需要建立三維立體模型來(lái)模擬污染氣體的擴(kuò)散過(guò)程?？紤]到三維模型的編程復(fù)雜度，而且污染氣體的擴(kuò)散在xoy平面上各向同性，可以將污染氣體在y方

21、向的擴(kuò)散等價(jià)為在x方向上的擴(kuò)散，此時(shí)便只需要建立xoz平面模型。如圖7所示。圖 7 污染氣體擴(kuò)散的等效二維示意圖對(duì)該模型的求解可以采用與問(wèn)題二同樣的分析方法，但要考慮到風(fēng)速是隨高度變化的，需要建立風(fēng)速與高度的關(guān)系式，并將其帶入到模型中?？紤]到風(fēng)速為0.6m/s ，1中給出的對(duì)應(yīng)的大氣穩(wěn)定度為E,對(duì)應(yīng)的風(fēng)速隨高度變化的關(guān)系式為： （9）其中，h 為高度， 為地面風(fēng)速。在本題中，將地面風(fēng)速看成0.6m/s .6.2 模型三的建立基于對(duì)問(wèn)題二的建模過(guò)程以及對(duì)問(wèn)題三的分析過(guò)程，現(xiàn)建立模型三如下： （10）其中，為x, y方向的擴(kuò)散系數(shù)， 為z方向的擴(kuò)散系數(shù)， 為時(shí)間，為污染源排放污染氣體的速度，為濃度

22、對(duì)時(shí)間的偏導(dǎo)，為濃度對(duì)x方向的偏導(dǎo)，為濃度對(duì)y方向的偏導(dǎo)，為對(duì)z方向的偏導(dǎo)，為濃度對(duì)x的二階偏導(dǎo)，為濃度對(duì)y的二階偏導(dǎo)。 分別表示風(fēng)速在x, y方向的分量。具體表達(dá)式為 （11）為了簡(jiǎn)化問(wèn)題，在求解過(guò)程中將z 方向的擴(kuò)散系數(shù)同樣當(dāng)作常數(shù)處理，即由于污染氣體在y方向的擴(kuò)散等價(jià)為在x方向上的擴(kuò)散，（10）式可以寫成 （12）邊界條件同樣取延拓邊界條件，即有 假設(shè)煙囪高度為10米，于是有源項(xiàng) 為由此得到在考慮風(fēng)向和高度影響時(shí)連續(xù)源的擴(kuò)散模型 （13）6.3 模型三的求解運(yùn)用ADI法求模型三的數(shù)值解，ADI法的格式為：根據(jù)（a）（b）的迭代格式及（13）中初邊值條件，運(yùn)用MATLAB 編程，可以得到任

23、意時(shí)刻的污染氣體濃度分布，其中截取60s、600s、1800s和3600s時(shí)濃度分布如圖8.圖 8（a）60s 時(shí)濃度分布圖圖 8（b）600s 時(shí)濃度分布圖圖 8（c）1800s 時(shí)濃度分布圖圖 8（d）3600s 時(shí)濃度分布圖圖 8 不同時(shí)刻濃度分布圖從圖8不同時(shí)刻濃度分布圖可以看出：任意時(shí)刻，煙囪口附近的污染氣體濃度最高；在x 方向，以煙囪口為中心，污染氣體的濃度在左右兩邊呈不對(duì)稱分布。開始時(shí)（t=60s），污染氣體在左邊擴(kuò)散到200米處，在右邊擴(kuò)散到580米處；隨著時(shí)間的推移，到t=600s時(shí)，污染氣體在左邊擴(kuò)散到100米處，在右邊仍為580米處；當(dāng)t =1800s時(shí)，污染氣體在左邊擴(kuò)

24、散到50米處，在右邊沒(méi)有什么變化，仍為580米處；當(dāng)t =3600s時(shí)，污染氣體在左邊擴(kuò)散至區(qū)域的邊緣處，右邊仍然在580米處。在z 方向，以煙囪口為中心，污染氣體的濃度在上下兩側(cè)同樣呈不對(duì)稱分布，任意時(shí)刻，下側(cè)的濃度擴(kuò)散范圍始終大于上側(cè)，到3600s 時(shí)，污染氣體已擴(kuò)散到距地面2m 左右。由此得出的結(jié)論是：由于風(fēng)向風(fēng)速及煙囪高度的影響，煙囪排放的污染氣體在空間中不是均勻分布的，在煙囪口處濃度最高，隨著時(shí)間的推移，上風(fēng)口處污染氣體的濃度及污染范圍變化不大，而在下風(fēng)口處，污染氣體的濃度和污染范圍在不斷地增大，到3600s 時(shí)，污染氣體已擴(kuò)散至近地面處?，F(xiàn)在將煙囪的高度設(shè)為20米，在其他條件不變的

25、情況下分析污染氣體的擴(kuò)散過(guò)程，并將其與煙囪高度為10米情況下進(jìn)行比較，截取3600s 時(shí)比較結(jié)果如下圖：圖 9（a）煙囪高度為10米的擴(kuò)散結(jié)果圖 9（b）煙囪高度為20米的擴(kuò)散結(jié)果圖 9 不同煙囪高度的擴(kuò)散結(jié)果從圖9的a ,b兩幅圖可以看出：當(dāng)煙囪高度為20米時(shí)，污染氣體在z 方向能向下擴(kuò)散5米左右，即地表及地表以上15米的范圍內(nèi)未受到污染氣體的影響。由此得出結(jié)論：增加煙囪高度能使污染氣體不擴(kuò)散至近地面，而是在十幾米的低空中因大氣運(yùn)動(dòng)擴(kuò)散出去，不會(huì)對(duì)地面人群造成影響。7、 給鋼鐵生產(chǎn)集團(tuán)治理污染氣體的建議報(bào)告目前，治理霧霾是人們最為關(guān)心的熱點(diǎn)問(wèn)題之一，霧霾形成的重要原因之一便是鋼鐵生產(chǎn)，為了減

26、小鋼鐵生產(chǎn)集團(tuán)煙囪污染氣體的排放對(duì)周邊地區(qū)大氣污染的影響，基于污染氣體擴(kuò)散模型，現(xiàn)提出以下建議：1. 減小污染氣體排放速度。當(dāng)排放量一定時(shí)，減小污染氣體排放速度，在單位時(shí)間內(nèi)大氣中污染物濃度的增量較小，減少了污染氣體在環(huán)境大氣中的積累，利于污染物的擴(kuò)散，使得大氣在較短時(shí)間內(nèi)恢復(fù)自凈能力。2. 考慮選址問(wèn)題。鋼鐵廠適宜建在城市的下風(fēng)口，因此在建廠之前需要調(diào)查了解城市的盛行風(fēng)向，將廠址選在遠(yuǎn)離居民區(qū)的下風(fēng)向，以此減小污染氣體對(duì)居民及周圍大氣的影響。3. 合理增加煙囪高度。在滿足國(guó)家有關(guān)煙囪高度標(biāo)準(zhǔn)的前提下，適當(dāng)增加煙囪高度，可以使得污染氣體的污染范圍保持在低空，不至于擴(kuò)散到近地面，從而保證了地面人

27、群不受污染氣體的影響。污染氣體的治理是一個(gè)長(zhǎng)期的過(guò)程，鋼鐵生產(chǎn)集團(tuán)在基于以上建議的基礎(chǔ)上，制訂出更好的治理管理措施，最終使得鋼鐵生產(chǎn)和環(huán)境治理達(dá)到雙贏的結(jié)果。8、 模型優(yōu)缺點(diǎn)分析8.1 模型優(yōu)點(diǎn)（1） 問(wèn)題一中建立了無(wú)風(fēng)及不考慮高度時(shí)的污染氣體擴(kuò)散模型。模型一建立過(guò)程中，考慮了污染氣體排放速度的實(shí)際情況，在二維擴(kuò)散方程中加入了有源項(xiàng)，對(duì)該二維擴(kuò)散方程的求解采用ADI法，得到迭代關(guān)系式，通過(guò)MATLAB編程，求出二維擴(kuò)散方程的數(shù)值解，并得到任意時(shí)刻污染氣體的濃度分布情況。通過(guò)附件一數(shù)據(jù)的檢驗(yàn)，該模型能夠很好的描述污染氣體的擴(kuò)散過(guò)程。（2） 問(wèn)題二在問(wèn)題一的基礎(chǔ)上考慮了風(fēng)向風(fēng)速的影響，在模型一的二

28、維擴(kuò)散方程中加入了風(fēng)向風(fēng)速的平流項(xiàng)，得到模型二。同樣采用ADI法，得到迭代關(guān)系式，利用MATLAB編程，求出模型二的數(shù)值解，并得到任意時(shí)刻污染氣體的濃度分布情況。通過(guò)附件二數(shù)據(jù)的檢驗(yàn)，該模型能夠很好的描述有風(fēng)情況下污染氣體的擴(kuò)散過(guò)程。（3） 問(wèn)題三在問(wèn)題二的基礎(chǔ)上考慮了高度的影響，在模型三的建立過(guò)程中，將復(fù)雜的三維空間擴(kuò)散模型轉(zhuǎn)化為二維擴(kuò)散模型，簡(jiǎn)化了編程的復(fù)雜度，根據(jù)模型三，利用MATLAB編程，模擬出了高度對(duì)污染氣體擴(kuò)散過(guò)程的影響。8.2 模型缺點(diǎn)（1） 在本文中將污染氣體擴(kuò)散系數(shù)當(dāng)作常數(shù)處理，未考慮擴(kuò)散系數(shù)與高度等因素的關(guān)系，導(dǎo)致模型與實(shí)際情況出現(xiàn)偏差。（2） 在本文中只考慮了同一類污染

29、氣體的擴(kuò)散過(guò)程，未考慮多種污染氣體同時(shí)存在時(shí)的擴(kuò)散情況，導(dǎo)致模型與實(shí)際出現(xiàn)偏差。9、 模型推廣與應(yīng)用模型一建立了無(wú)風(fēng)及不考慮高度時(shí)的污染氣體擴(kuò)散模型，該模型考慮了污染氣體排放速度的實(shí)際情況，在二維擴(kuò)散方程中加入了有源項(xiàng)，采用ADI法得到擴(kuò)散方程的迭代式，利用MATLAB編程，求出其數(shù)值解，得到任意時(shí)刻污染氣體的濃度分布。經(jīng)檢驗(yàn)該模型與實(shí)際相符，所以該模型可用于求解平面上無(wú)風(fēng)情況下的污染氣體擴(kuò)散過(guò)程。模型二建立了有風(fēng)情況下的污染氣體擴(kuò)散模型，采用與模型一同樣的處理方法，得到任意時(shí)刻污染氣體的濃度分布。經(jīng)檢驗(yàn)該模型與實(shí)際相符，所以該模型可用于求解平面上有風(fēng)情況下的污染氣體擴(kuò)散過(guò)程。模型三綜合考慮了

30、風(fēng)向風(fēng)速及高度等因素的影響，使模型能夠更好的運(yùn)用于實(shí)際情況中，能較好的模擬真實(shí)情況下污染氣體的擴(kuò)散過(guò)程，為鋼鐵生產(chǎn)集團(tuán)提供理論依據(jù)。參考文獻(xiàn)1 王志春 宋麗莉 何秋生 劉愛(ài)君 劉榮 葉燕翔，風(fēng)速隨高度變化的曲線模型分析，熱帶氣象學(xué)報(bào)，第23卷第6期：690-692頁(yè)，20072 姜啟源 謝金星 葉俊，數(shù)學(xué)模型（第四版），北京：高等教育出版社，20113 Rafael B. Storch Luiz C.G. Pimentel Helcio R.B. Orlande，Identification of atmospheric boundary layer parameters by inverse

31、 problem , atmospheric environment 41 :1417-1425,2007附錄：1. 無(wú)風(fēng)模型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 intervalsN=60;xb=0.0;xe=1000.0;yb=0.0;ye=1000.0;J=101;K=101;% for time intervalst0=0;dt=1;tf=dt*N;% the number of mesh po

32、intsJ01=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/(dx2);nuy=dt*a/(dy2);%to obtain the coordinatesx=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; endend%-Big loop for time t -for n =1:N t1 = t + dt;t2 = t + d

33、t/2;%-sweep in x-direction -%boundary conditionsu2(:,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

34、(1)=0;%these are for zero boundary conditions b(J)=0; ut = Ab;%solve the tri-diagonal matrix. u2(:,k)=ut;end %finishi x-sweep%-sweep in y-direction -%boundary conditionsu1(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)+

35、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 = Ab;%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(

36、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(原始數(shù)據(jù)曲線圖)suptitle(t=60s)2. 無(wú)風(fēng)時(shí)源項(xiàng)ffunction result = f(x,y)if (x-500)2+(y-500)2=1 result=300;else result=0;endend3. 有風(fēng)模型clc;clear;% N is the

37、total time steps% J is the total mesh points in X% K is the total mesh points in Y % for space intervalsN=360;xb=0.0;xe=1000.0; yb=0.0;ye=1000.0;J=101;K=101;% for time intervalst0=0;dt=1;tf=dt*N; % the number of mesh pointsJ01=J-1;K01=K-1; dx=(xe-xb)/(J-1);dy=(ye-yb)/(K-1);%dt=(tf-t0)/(N-1); a=10;m1

38、=sqrt(2)*0.3;m2=sqrt(2)*0.3;nux=dt*a/(dx2);nuy=dt*a/(dy2);mux=dt*m1/(2*dx);muy=dt*m2/(2*dy);%to obtain the coordinatesx=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; endend%-Big loop for time t -for n =1:N t1 = t + dt;t2 = t + dt/2;%-sweep in x

39、-direction -%boundary conditionsu2(:,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; en

40、d A(1,1)=1; A(J,J)=1; b(1)=0;%these are for zero boundary conditions b(J)=0; ut = Ab;%solve the tri-diagonal matrix. u2(:,k)=ut;end %finishi x-sweep%-sweep in y-direction -%boundary conditionsu1(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)

41、= (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 = Ab;%solve the tri-diagonal matrix. u1(j,:)=ut; end %finish

42、i x-sweep endmesh(x,y,u1)4. 有風(fēng)時(shí)的源項(xiàng)f2function result = f2(x,z)if (x-500)2+(z-10)2=1 result=300;else result=0;endEnd5. 原始數(shù)據(jù)濃度的導(dǎo)出clc;clear;load (a180.txt);c1=a180(:,3);dlmwrite(a1801s.txt,c1)6. 考慮風(fēng)速風(fēng)向和高度的模型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 intervalsN=60;xb=0.0;xe=1000.0;zb=0.0;ze=100.0;J=101;K=101;% for time intervalst0=0;dt=1;tf=dt*N;% the number of mesh pointsJ01=J-