空氣污染數(shù)學建模.doc

A.污染氣體的傳播擴散摘要鋼鐵生產(chǎn)排放的污染氣體是造成霧霾的重要原因之一，研究污染氣體的擴散特征，正確模擬污染氣體的擴散過程，能夠為鋼鐵生產(chǎn)集團提出更好的治理管理措施，具有實際意義。針對問題一：污染氣體的排放速度為300m/s,在不考慮風向風速及高度影響的情況下，此問題即為二維平面的連續(xù)點源擴散問題，由此在二維xoy 平面上建立連續(xù)點源擴散方程模型，其中 為氣體擴散系數(shù)，本文中取為常數(shù)10，f(x,y,t ) 為污染氣體的排放速度，在本文中恒為300m/s ；對上述偏微分方程模型，本文采用ADI法（Alternating direction implicit，交替方向隱式法）求解出迭代格式，利用MATLAB編程，求出模型一的數(shù)值解，并得到任意時刻污染氣體的濃度分布情況。通過SPSS軟件，對附件一所給的原始實際數(shù)據(jù)與模型一求解得到的模擬值進行顯著性檢驗，檢驗結(jié)果顯示該模型與實際情況吻合。針對問題二：考慮風向風速對污染氣體擴散過程的影響時，在基于對問題一求解的基礎(chǔ)上，在模型一的擴散方程模型中加入風向風速的平流項，由此得到有風情況下的模型，其中 分別為風速在x, y方向的分量；對此模型同樣采用ADI法求出迭代格式，利用MATLAB編程，求出模型二的數(shù)值解，并得到任意時刻污染氣體的濃度分布情況。通過SPSS軟件，對附件二所給的原始實際數(shù)據(jù)與模型二求解得到的模擬值進行顯著性檢驗，檢驗結(jié)果顯示該模型與實際情況吻合。針對問題三：考慮有風時增加高度的影響，此問題即為三維空間的污染氣體擴散問題，考慮到三維模型的編程復雜度，而且污染氣體的擴散在xoy平面上各向同性，可以將污染氣體在y方向的擴散等價為在x方向上的擴散，此時便只需要建立xoz平面上的擴散模型。在基于對問題二求解的基礎(chǔ)上，在模型二的擴散方程中增加高度項，由此得到模型三為，其中 為z 方向的擴散系數(shù)；對該擴散方程同樣采用ADI法求出迭代格式，利用MATLAB編程，求出模型二的數(shù)值解，并得到任意時刻污染氣體的濃度分布情況。關(guān)鍵詞：污染氣體 擴散方程 ADI法 數(shù)值解1、 問題重述目前，治理霧霾是人們最為關(guān)心的熱點問題之一。中國社科院發(fā)布的氣候變化綠皮書中提及，霧霾形成的原因里，重工業(yè)、車輛尾氣、土方施工都榜上有名，其中鋼鐵生產(chǎn)也是造成霧霾的重要原因之一。某鋼鐵生產(chǎn)集團煙囪污染氣體的排放對周邊地區(qū)大氣污染的影響非常大，為了提出更好的治理管理措施，需要對其污染氣體擴散的特征進行分析。現(xiàn)在，我們需要在三種情況下考慮污染氣體的擴散過程：1. 在不考慮風向和高度影響的情況下，建立模型，模擬某鋼鐵生產(chǎn)集團的煙囪排放污染氣體的擴散過程，假設(shè)煙囪的排放速度為300m/s。2. 考慮風向為東北風，平均風速0.6m/s的情況下，模擬污染氣體的傳播擴散過程。3. 在考慮風向的基礎(chǔ)上增加高度的影響，建立模型，模擬污染氣體的傳播擴散過程。4. 基于上述模型結(jié)論，給該鋼鐵生產(chǎn)集團提供一個污染氣體治理建議報告。2、 問題分析鋼鐵生產(chǎn)集團煙囪污染氣體的排放對周邊地區(qū)大氣污染的影響非常大，為了提出更好的治理管理措施，需要對其污染氣體擴散的特征進行分析。污染氣體擴散是與時間、空間相關(guān)的連續(xù)性問題，本文利用微分方程建立模型，通過ADI法（Alternating direction implicit，交替方向隱式法）求解微分方程的迭代式，利用MATLAB編程求其數(shù)值解來模擬擴散過程。對于問題一，在不考慮風速和高度影響的情況下，污染氣體的擴散問題即為二維平面點源的擴散問題；考慮到煙囪的排放速度為300m/s，需要在二維擴散方程中加入有源項，由此建立擴散模型。為檢驗模型的準確性，需利用SPSS軟件，將擴散方程的數(shù)值解與實際數(shù)據(jù)進行顯著性分析。對于問題二，在問題一的基礎(chǔ)上，考慮了風向風速的影響，需要在模型一中加入風向風速對應(yīng)的平流項，由此建立模型二。為檢驗模型的準確性，同樣需要利用SPSS軟件，將模型二的數(shù)值解與實際數(shù)據(jù)進行顯著性分析。對于問題三，基于對問題一、二的建模思路，首先建立不考慮風速風向只考慮高度時的模型，此時應(yīng)分析擴散系數(shù)與高度的關(guān)系，建立模型三；然后在模型三的基礎(chǔ)上考慮風速風向的影響，此時應(yīng)分析風速與高度的關(guān)系，建立模型四。為使模型更能符合實際情況下污染氣體的擴散過程，本問題還考慮了大氣靜力穩(wěn)定度對擴散過程的影響。綜上所述，本問題可以看成是求解偏微分方程中擴散方程的數(shù)值解問題。3、 模型假設(shè)與符號說明1. 假設(shè)煙囪的排放速度恒為300m/s；2. 假設(shè)污染氣體中只含有同一類污染物，污染氣體的擴散系數(shù)相同；3. 假設(shè)不考慮氣體受到的重力和浮力；4. 假設(shè)在整個擴散過程中污染氣體不發(fā)生沉降、分解，不發(fā)生化學反應(yīng)；5. 假設(shè)地面以及地標地物對氣體無吸收；6. 假設(shè)在有風情況下，風向水平，風速風向恒定；u : 污染物濃度f : 污染物排放速度v : 風速 : 水平方向擴散系數(shù)： 風速在x 方向的分量： 風速在y 方向的分量 ： 垂直方向擴散系數(shù)4、 對問題一的分析與建模4.1 問題分析問題一在不考慮風速和高度影響下，煙囪以恒定速度排放污染氣體，在無限大平面上，將煙囪口看成是點源，此問題即為二維平面點源的擴散問題，建立二維擴散方程模型，通過ADI法求其數(shù)值解，即可得到任意時刻污染氣體的濃度分布。4.2 模型一的建立本文將煙囪口當作二維平面上的點源來考慮，分析附件allresults_1所給的數(shù)據(jù)，可以得出的結(jié)論有：污染氣體擴散的平面圖為10001000（單位：m）的正方形區(qū)域；以點（500,500）為中心，半徑為1m的領(lǐng)域內(nèi)污染氣體的濃度最高，據(jù)此確定煙囪口的位置坐標為（500,500）。假設(shè)污染源為連續(xù)點源，建立有源項的二維擴散方程： （1）其中，為擴散系數(shù)， 為時間，為污染源排放污染氣體的速度，為濃度對時間的偏導， 為濃度對x的二階偏導，為濃度對y的二階偏導。現(xiàn)在考慮初邊值條件，初始時刻，平面上各點處濃度均為0，即有邊界條件取延拓邊界條件，即有 根據(jù)題意知，在污染源處，污染物排放速度為300m/s，假設(shè)煙囪口是以點（500,500）為中心，1m 為半徑的圓形區(qū)域，在10001000的平面內(nèi)，上述圓形區(qū)域任然可以當作點源來處理，在圓形區(qū)域內(nèi)有源項為300，在圓形區(qū)域外有源項為0，即有由此得到在不考慮風向和高度影響時連續(xù)源的擴散模型 （2）另外，根據(jù)中國大百科全書大氣擴散卷，選定擴散系數(shù)。4.3 模型一求解本文采用ADI法求解（1）式數(shù)值解，將（1）式轉(zhuǎn)化為ADI格式為 （a）（b）對（a）式進行轉(zhuǎn)換，令： 并將帶入（a）式中，得： (3)將(3)式寫成矩陣方程組形式： 同樣對（b）式進行轉(zhuǎn)換，令，并將帶入（a）式中，得： (4)將(4)式寫成矩陣方程組形式： 根據(jù)（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個單位，遠遠高于周圍其他區(qū)域，整個濃度分布區(qū)域近似為一條垂直于xoy平面的直線；隨著時間的推移，到t=600s時，煙囪口處污染氣體的濃度近似為1100個單位，且周圍其他區(qū)域內(nèi)的濃度有所增大，整個濃度分布區(qū)域由60s的直線狀變?yōu)槌噬霞庀聢A的塔狀；當t =1800s時，煙囪口處污染氣體的濃度近似為1300個單位；當t =3600s時，煙囪口處污染氣體的濃度已達到1700個單位左右，從“塔”的形狀來看，“塔尖”在不斷地增高，“塔底”在不斷地向外擴展，到3600s 時已擴散到半徑為800米的圓形區(qū)域內(nèi)。4.4 模型一的檢驗首先通過MATLAB畫圖，比較模型一的濃度分布圖與用實際數(shù)據(jù)畫出的濃度分布圖的區(qū)別，截取120s 和3600s 的濃度對比圖如下：圖 2（a）120s 時模擬值與實際數(shù)據(jù)濃度分布對比圖圖 2（b）3600s 時模擬值與實際數(shù)據(jù)濃度分布對比圖圖 2 不同時刻模擬值與實際數(shù)據(jù)濃度分布對比圖由圖2（a）（b）兩幅圖可以看出：由模型一得到的濃度分布圖與用實際數(shù)據(jù)畫出的濃度分布圖在濃度大小和濃度分布形狀上近似一致。為了更進一步的分析兩者之間的誤差，現(xiàn)在運用SPSS進一步做模型檢驗。取120s 和3600s 時由模型一得到的濃度分布數(shù)據(jù)與實際數(shù)據(jù)進行顯著性檢驗。首先畫出模擬值與實際數(shù)據(jù)的比較圖如下：圖 3（a）120s 模擬值與實際數(shù)據(jù)比較圖圖 3（b）3600s 時模擬值與實際數(shù)據(jù)比較圖圖 3 不同時刻模擬值與實際數(shù)據(jù)比較圖從圖3可以直觀的看出，模擬值與實際數(shù)據(jù)之間成明顯的線性關(guān)系，現(xiàn)在進一步對模擬值與實際數(shù)據(jù)做線性回歸分析，取120s時的模擬數(shù)據(jù)和實際數(shù)據(jù)做線性回歸分析，SPSS運行結(jié)果如下：Anovaa模型平方和df均方FSig.1回歸12176482.538112176482.5382966236.658.000b殘差41867.174101994.105總計12218349.71110200a. 因變量: V1b. 預測變量: (常量), V2。系數(shù)a模型非標準化系數(shù)標準系數(shù)tSig.B 的 95.0% 置信區(qū)間B標準 誤差試用版下限上限實際值(常量)-.089.020-4.432.000-.129-.050模擬值1.260.001.9981722.277.0001.2591.262a. 因變量: V1從運行結(jié)果中可以看出：回歸檢驗的p 值為0.000，小于0.05，拒絕原假設(shè)，即模擬值與實際數(shù)據(jù)線性回歸顯著。且實際值與模擬值之比為1.260，由此可得出模擬值與實際數(shù)據(jù)基本吻合，即可以用模型一來模擬無風條件下污染氣體的擴散過程。5、 對問題二的分析與建模5.1 問題分析問題二在問題一的基礎(chǔ)上增加了風速和風向的影響，需要在模型一中加入風向風速對應(yīng)的平流項，同樣通過ADI法求其數(shù)值解，得到任意時刻污染氣體的濃度分布。5.2 模型二的建立假設(shè)只考慮水平風向，且風向恒為東北風，風速恒為0.6m/s，將東北風向分解到x ,y 方向上，有其中， 分別為x ,y 方向上風速分量。在模型一的基礎(chǔ)上加入風速風向?qū)?yīng)的平流項，即有 （5）其中，為擴散系數(shù)， 為時間，為污染源排放污染氣體的速度，為濃度對時間的偏導，為濃度對x的偏導，為濃度對y的偏導，為濃度對x的二階偏導，為濃度對y的二階偏導。邊界條件取延拓邊界條件，即有 有源項 與問題一中相同，即由此得到在考慮風向影響時連續(xù)源的擴散模型 （6）5.3 模型二的求解運用ADI法求模型二的數(shù)值解，ADI法的格式為：令 ，并將帶入（a）式中,對（a）式進行變形，得：（7）將（7）式轉(zhuǎn)化為矩陣方程組形式 k=1,2,3,.,K令，并將帶入（a）式中,對（b）式進行變形，得：（8）同樣，可以將（8）式轉(zhuǎn)化為矩陣方程組形式 J=1,2,3,.J根據(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個單位，遠遠高于周圍其他區(qū)域，整個濃度分布區(qū)域近似為一條垂直于xoy平面的直線；隨著時間的推移，到t=600s時，煙囪口處污染氣體的濃度近似為900個單位，且在 （即下風口）區(qū)域內(nèi)的濃度有所增大，整個濃度分布區(qū)域由60s的直線狀變?yōu)榘脲F形；當t =1800s和t =3600s時，煙囪口處污染氣體的濃度任然近似為900個單位，但半錐形底部半圓的面積在不斷擴大，到3600s 時已擴散到整個的區(qū)域內(nèi)。由此得出的結(jié)論是：由于風向風速的存在，煙囪排放的污染氣體在平面中不是均勻分布的，在煙囪口處濃度最高，且煙囪口處污染氣體濃度隨時間增加不明顯，在下風口污染氣體濃度隨時間增加而增大，且污染范圍隨時間增加而擴大。5.4 模型二的檢驗首先通過MATLAB畫圖，比較模型一的濃度分布圖與用實際數(shù)據(jù)畫出的濃度分布圖的區(qū)別，截取120s 和3600s 的濃度對比圖如下：圖 5（a）120s 時模擬值與實際數(shù)據(jù)濃度分布對比圖圖 5（b）3600s 時模擬值與實際數(shù)據(jù)濃度分布對比圖圖 5 不同時刻模擬值與實際數(shù)據(jù)濃度分布對比圖由圖5（a）（b）兩幅圖可以看出：由模型二得到的濃度分布圖與用實際數(shù)據(jù)畫出的濃度分布圖在濃度大小和濃度分布形狀上近似一致。為了更進一步的分析兩者之間的誤差，現(xiàn)在運用SPSS進一步做模型檢驗。取120s 和3600s 時由模型二得到的濃度分布數(shù)據(jù)與實際數(shù)據(jù)進行顯著性檢驗。首先畫出模擬值與實際數(shù)據(jù)的比較圖如下：圖 6（a）120s 模擬值與實際數(shù)據(jù)比較圖圖 6（b）3600s 模擬值與實際數(shù)據(jù)比較圖圖 6 不同時刻模擬值與實際數(shù)據(jù)比較圖從圖6可以直觀的看出，模擬值與實際數(shù)據(jù)之間成明顯的線性關(guān)系，現(xiàn)在進一步對模擬值與實際數(shù)據(jù)做線性回歸分析，取120s時的模擬數(shù)據(jù)和實際數(shù)據(jù)做線性回歸分析，SPSS運行結(jié)果如下：Anovaa模型平方和df均方FSig.1回歸5789448.83515789448.8351428052.592.000b殘差41347.629101994.054總計5830796.46410200a. 因變量: V1b. 預測變量: (常量), V2。系數(shù)a模型非標準化系數(shù)標準系數(shù)tSig.B標準 誤差試用版實際值(常量)-.074.020-3.715.000模擬值1.281.001.9961195.012.000a. 因變量: V1從運行結(jié)果中可以看出：回歸檢驗的p 值為0.000，小于0.05，拒絕原假設(shè)，即模擬值與實際數(shù)據(jù)線性回歸顯著。且實際值與模擬值之比為1.281，由此可得出模擬值與實際數(shù)據(jù)基本吻合，即可以用模型二來模擬有風條件下污染氣體的擴散過程。6、 對問題三的分析與建模6.1 問題分析問題三在問題二的基礎(chǔ)上增加了高度的影響，假設(shè)煙囪高度為10米，此時二維平面模型已不再適用，需要建立三維立體模型來模擬污染氣體的擴散過程?？紤]到三維模型的編程復雜度，而且污染氣體的擴散在xoy平面上各向同性，可以將污染氣體在y方向的擴散等價為在x方向上的擴散，此時便只需要建立xoz平面模型。如圖7所示。圖 7 污染氣體擴散的等效二維示意圖對該模型的求解可以采用與問題二同樣的分析方法，但要考慮到風速是隨高度變化的，需要建立風速與高度的關(guān)系式，并將其帶入到模型中?？紤]到風速為0.6m/s ，1中給出的對應(yīng)的大氣穩(wěn)定度為E,對應(yīng)的風速隨高度變化的關(guān)系式為： （9）其中，h 為高度， 為地面風速。在本題中，將地面風速看成0.6m/s .6.2 模型三的建立基于對問題二的建模過程以及對問題三的分析過程，現(xiàn)建立模型三如下： （10）其中，為x, y方向的擴散系數(shù)， 為z方向的擴散系數(shù)， 為時間，為污染源排放污染氣體的速度，為濃度對時間的偏導，為濃度對x方向的偏導，為濃度對y方向的偏導，為對z方向的偏導，為濃度對x的二階偏導，為濃度對y的二階偏導。 分別表示風速在x, y方向的分量。具體表達式為 （11）為了簡化問題，在求解過程中將z 方向的擴散系數(shù)同樣當作常數(shù)處理，即由于污染氣體在y方向的擴散等價為在x方向上的擴散，（10）式可以寫成 （12）邊界條件同樣取延拓邊界條件，即有 假設(shè)煙囪高度為10米，于是有源項 為由此得到在考慮風向和高度影響時連續(xù)源的擴散模型 （13）6.3 模型三的求解運用ADI法求模型三的數(shù)值解，ADI法的格式為：根據(jù)（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時，污染氣體在左邊擴散至區(qū)域的邊緣處，右邊仍然在580米處。在z 方向，以煙囪口為中心，污染氣體的濃度在上下兩側(cè)同樣呈不對稱分布，任意時刻，下側(cè)的濃度擴散范圍始終大于上側(cè)，到3600s 時，污染氣體已擴散到距地面2m 左右。由此得出的結(jié)論是：由于風向風速及煙囪高度的影響，煙囪排放的污染氣體在空間中不是均勻分布的，在煙囪口處濃度最高，隨著時間的推移，上風口處污染氣體的濃度及污染范圍變化不大，而在下風口處，污染氣體的濃度和污染范圍在不斷地增大，到3600s 時，污染氣體已擴散至近地面處?，F(xiàn)在將煙囪的高度設(shè)為20米，在其他條件不變的情況下分析污染氣體的擴散過程，并將其與煙囪高度為10米情況下進行比較，截取3600s 時比較結(jié)果如下圖：圖 9（a）煙囪高度為10米的擴散結(jié)果圖 9（b）煙囪高度為20米的擴散結(jié)果圖 9 不同煙囪高度的擴散結(jié)果從圖9的a ,b兩幅圖可以看出：當煙囪高度為20米時，污染氣體在z 方向能向下擴散5米左右，即地表及地表以上15米的范圍內(nèi)未受到污染氣體的影響。由此得出結(jié)論：增加煙囪高度能使污染氣體不擴散至近地面，而是在十幾米的低空中因大氣運動擴散出去，不會對地面人群造成影響。7、 給鋼鐵生產(chǎn)集團治理污染氣體的建議報告目前，治理霧霾是人們最為關(guān)心的熱點問題之一，霧霾形成的重要原因之一便是鋼鐵生產(chǎn)，為了減小鋼鐵生產(chǎn)集團煙囪污染氣體的排放對周邊地區(qū)大氣污染的影響，基于污染氣體擴散模型，現(xiàn)提出以下建議：1. 減小污染氣體排放速度。當排放量一定時，減小污染氣體排放速度，在單位時間內(nèi)大氣中污染物濃度的增量較小，減少了污染氣體在環(huán)境大氣中的積累，利于污染物的擴散，使得大氣在較短時間內(nèi)恢復自凈能力。2. 考慮選址問題。鋼鐵廠適宜建在城市的下風口，因此在建廠之前需要調(diào)查了解城市的盛行風向，將廠址選在遠離居民區(qū)的下風向，以此減小污染氣體對居民及周圍大氣的影響。3. 合理增加煙囪高度。在滿足國家有關(guān)煙囪高度標準的前提下，適當增加煙囪高度，可以使得污染氣體的污染范圍保持在低空，不至于擴散到近地面，從而保證了地面人群不受污染氣體的影響。污染氣體的治理是一個長期的過程，鋼鐵生產(chǎn)集團在基于以上建議的基礎(chǔ)上，制訂出更好的治理管理措施，最終使得鋼鐵生產(chǎn)和環(huán)境治理達到雙贏的結(jié)果。8、 模型優(yōu)缺點分析8.1 模型優(yōu)點（1） 問題一中建立了無風及不考慮高度時的污染氣體擴散模型。模型一建立過程中，考慮了污染氣體排放速度的實際情況，在二維擴散方程中加入了有源項，對該二維擴散方程的求解采用ADI法，得到迭代關(guān)系式，通過MATLAB編程，求出二維擴散方程的數(shù)值解，并得到任意時刻污染氣體的濃度分布情況。通過附件一數(shù)據(jù)的檢驗，該模型能夠很好的描述污染氣體的擴散過程。（2） 問題二在問題一的基礎(chǔ)上考慮了風向風速的影響，在模型一的二維擴散方程中加入了風向風速的平流項，得到模型二。同樣采用ADI法，得到迭代關(guān)系式，利用MATLAB編程，求出模型二的數(shù)值解，并得到任意時刻污染氣體的濃度分布情況。通過附件二數(shù)據(jù)的檢驗，該模型能夠很好的描述有風情況下污染氣體的擴散過程。（3） 問題三在問題二的基礎(chǔ)上考慮了高度的影響，在模型三的建立過程中，將復雜的三維空間擴散模型轉(zhuǎn)化為二維擴散模型，簡化了編程的復雜度，根據(jù)模型三，利用MATLAB編程，模擬出了高度對污染氣體擴散過程的影響。8.2 模型缺點（1） 在本文中將污染氣體擴散系數(shù)當作常數(shù)處理，未考慮擴散系數(shù)與高度等因素的關(guān)系，導致模型與實際情況出現(xiàn)偏差。（2） 在本文中只考慮了同一類污染氣體的擴散過程，未考慮多種污染氣體同時存在時的擴散情況，導致模型與實際出現(xiàn)偏差。9、 模型推廣與應(yīng)用模型一建立了無風及不考慮高度時的污染氣體擴散模型，該模型考慮了污染氣體排放速度的實際情況，在二維擴散方程中加入了有源項，采用ADI法得到擴散方程的迭代式，利用MATLAB編程，求出其數(shù)值解，得到任意時刻污染氣體的濃度分布。經(jīng)檢驗該模型與實際相符，所以該模型可用于求解平面上無風情況下的污染氣體擴散過程。模型二建立了有風情況下的污染氣體擴散模型，采用與模型一同樣的處理方法，得到任意時刻污染氣體的濃度分布。經(jīng)檢驗該模型與實際相符，所以該模型可用于求解平面上有風情況下的污染氣體擴散過程。模型三綜合考慮了風向風速及高度等因素的影響，使模型能夠更好的運用于實際情況中，能較好的模擬真實情況下污染氣體的擴散過程，為鋼鐵生產(chǎn)集團提供理論依據(jù)。參考文獻1 王志春 宋麗莉 何秋生 劉愛君 劉榮 葉燕翔，風速隨高度變化的曲線模型分析，熱帶氣象學報，第23卷第6期：690-692頁，20072 姜啟源 謝金星 葉俊，數(shù)學模型（第四版），北京：高等教育出版社，20113 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 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 pointsJ01=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 + dt/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(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)+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(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. 無風時源項ffunction result = f(x,y)if (x-500)2+(y-500)2=1 result=300;else result=0;endend3. 有風模型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=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=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-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; end 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) = (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 %finishi x-sweep endmesh(x,y,u1)4. 有風時的源項f2function result = f2(x,z)if (x-500)2+(z-10)2=1 result=300;else result=0;endEnd5. 原始數(shù)據(jù)濃度的導出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 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-