基于多相質點網格法的大氣污染仿真模擬

于瀟萌， 曹 樂*， 嚴家德， 王巍巍

(1.南京信息工程大學中國氣象局氣溶膠與云降水重點開放實驗室，南京 210044；2.南京信息工程大學中國 氣象局綜合氣象觀測實習基地，南京 210044；3.南京信息工程大學大氣與環境實驗教學中心，南京 210044)

近年來，隨著城市化與工業化進程的不斷推進，大氣污染問題日益嚴重。為了更好地完成空氣污染問題的綜合治理，學者們在大氣環境領域做了諸多研究[1]，以預測空氣污染物的散布及變化規律。這些研究包括對下墊面條件的處理、大氣擴散與大氣輸送的尺度的拓展、多層次污染之間的相互作用等多個方面，其研究成果也被廣泛應用于城市規劃、生態建設、信息通信、醫療健康等不同領域。而當前大氣環境污染相關研究工作的瓶頸之一，在于如何準確模擬大氣污染物的擴散規律[2]。大氣污染物擴散規律制約眾多，包括了多變的氣象條件、污染源排放特性、下墊面邊界條件、復雜物理化學過程等因素，因此學者們開發出了眾多有著不同適用范圍的擴散模型，以應對復雜的影響因素。而選取適當的擴散模型有助于更好地提高擴散模擬的質量，改進區域環境質量的預報效果[3-6]。

目前應用相對廣泛的大氣擴散模型大致可分為兩類[7]，一類是基于高斯理論的半經驗模型，另一類是計算流體動力學(computational fluid dynamics，CFD)模型。相比較簡單的高斯經典擴散模式，CFD模型可以更好地刻畫湍流以及復雜地形條件對污染物擴散造成的影響[8]，因此，CFD模擬成為了大氣污染模擬的主要手段[9]，并被廣泛應用于不同尺度的大氣流動和擴散問題上[10]。Rahimi等[11]用CFD模型模擬了伊朗伊斯法罕煉油廠不同煙囪釋放的氣體污染物的擴散，比較并考慮了3種不同湍流模型(k-ε模型、RNGk-ε模型、realizablek-ε模型)對擴散造成的影響，并將模擬結果與煉油廠實際觀測數據進行對比驗證。Jeanjean等[12]利用OpenFOAM軟件平臺，采用k-ε模型對大氣污染物濃度進行了CFD模擬，以研究樹木與氣流動力學之間的相互作用。Tominaga等[13]通過比較不同羽流浮力條件下的RANS模擬結果，與實驗結果對比驗證，探究使用Boussinesq近似進行RANS(reynolds-averaged navier-stokes)模擬的可行性。楊光俊等[14]為研究燃煤電廠的煙氣擴散，采用CFD方法，對燃煤電廠兩種排放方式(煙塔合一和煙囪結構)的擴散形態進行模擬。

然而，當前諸多應用CFD軟件的數值模擬研究，多數聚焦于不同排放源諸如火電廠、冷卻塔、汽車尾氣等對污染擴散造成的影響，或關注于空氣污染物輸送過程中化學轉化、干沉降和濕沉降清除等對污染擴散的作用，對于湍流作用下污染物擴散具體形態變化及濃度變化定量分析的研究還比較少。而離散粒子法(discrete particle model，DPM)[15]經過學者們的驗證，可以準確模擬湍流作用下顆粒污染物的擴散形態變化與濃度變化情況。榮雙等[16]應用DPM方法研究街道峽谷內不同高度處污染物濃度的影響，結果表明風向與街谷長寬比的變化都會對污染物的濃度產生影響。付志民等[17]為研究街道峽谷的不同幾何結構對風場與污染物的影響，采用DPM方法進行模擬，結果表明遞增型街谷模型更有利于顆粒污染物的擴散。

但是，傳統DPM方法在計算粒子之間碰撞作用的計算量極大，難以應用到范圍較大的污染物擴散模擬中，因此只適用于室內或短距離之間的污染物擴散模擬。而以DPM為基礎發展起來的多相質點網格方法(multi-phase particle-in-cell，MP-PIC)[18]以固相應力的方式取代顆粒之間的碰撞過程，可以大大提高計算效率，也為提高拉格朗日粒子模型的預報尺度提供了可能。因此，使用以MP-PIC為基礎的拉格朗日粒子追蹤模型，對理想條件下大氣污染擴散問題進行大渦仿真模擬，從而分析污染物擴散形態變化特征及污染物數濃度時空分布變化特征，以探究不同風速條件下湍流對污染物擴散造成的影響。這一工作無論在大氣環境模式發展，還是空氣質量預報業務，亦或者關于大氣污染潛在影響的研究中都具有相當的應用價值。

1 研究方法與算例設置

本文所使用的計算平臺為開源軟件OpenFOAM(open field operation and manipulation)[19]，并利用基于MP-PIC方法的拉格朗日粒子模型對單點源污染物不同風速條件下的擴散情況進行模擬計算。下面對軟件平臺及具體求解器進行詳細介紹。

1.1 求解器介紹

OpenFOAM作為一款基于Linux 平臺的開源軟件，其本質是一個功能強大的預編譯功能庫。它自帶許多開發完備的基礎算例與現成的求解器，用戶可以通過命令直接調用求解器對具體問題進行求解，或根據需要對這些功能庫進行修改以滿足自己的需求[20-21]。這些標準求解器包括了不可壓縮流求解器、可壓縮流求解器、多相流求解器等。而本研究選用標準求解器MPPICFoam作為求解器，該求解器為瞬態流體-固體粒子求解器，可以將粒子打包為粒子云，同時附帶離散相對連續相的影響，適合研究氣固兩相流之間的相互作用。

MPPICFoam中所使用的MP-PIC方法于1996年被提出，最初被應用于模擬密相顆粒流的工作。該方法將流體視為連續相、顆粒視為離散相進行求解，分別采用歐拉坐標系與拉格朗日坐標系進行描述，以追蹤粒子的具體軌跡[22-23]。與傳統DPM方法相比，MP-PIC方法用固相應力的方式取代顆粒之間的碰撞過程，可以大大提高計算效率。

在MP-PIC方法中，氣相連續性方程為

(1)

式(1)中：θg表示氣相體積分數；ρg表示氣相密度，kg/m3；ug表示氣相速度，m/s。

而氣相動量方程可表示為

·θgτg+θgρgg-F

(2)

式(2)中：τg表示氣相應力張量，Pa；p表示壓強，Pa；F表示氣固相間動量交換力，N/m3，即顆粒所受到的拖曳力；g表示重力加速度，m/s2。τg的表達式為

(3)

式(3)中：μ1為分子黏度，kg/(m·s)；μt為湍流黏度，kg/(m·s)，可以由次網格模型計算得出；S表示形變率張量，s-1；I為單位張量；S可表示為

(4)

相比于氣相方程，固相粒子在空間中的分布情況可表示為[24]

(5)

式(5)中：f表示粒子分布函數；up表示了平均粒子速度，m/s；A表示粒子加速度，m/s2，其表達式為

(6)

式(6)中：Dp表示阻力系數，s-1，θs、ρs分別表示固相體積分數、固相平均密度，kg/m3；ps表示固相應力，Pa。其中，固相體積分數θs與氣相體積分數θg之間滿足關系式：

θg+θs=1

(7)

而固相應力ps的表達式為

(8)

氣相與固相之間通過氣固相間動量交換力F(單位為N/m3)進行耦合，F與粒子分布函數f之間的關系為

F=?fVp[βp(ug-up)-p]dVpdρsdup

(9)

式(9)中：Vp為氣相體積，m3；βp表示曳力系數，kg/(m3·s)，具體公式可參考文獻[24]。

本次模擬中湍流模型采用大渦模擬(large eddy simulation，LES)，應用的次網格模型為Smagorinsky模型[25]，次網格尺度渦黏計算公式為

(10)

式(10)中：Cs表示模型常數；Δ為過濾尺寸長度，

表示過濾后的應變率，s-1，表達式為

(11)

1.2 算例基本設置

圖1所示為本次模擬的計算區域設置，整個計算區域形狀為長方體，大小為1 500 m×200 m×200 m。左側為計算域入口，而污染源距離左側入口約100 m，污染源高度為10 m，對外不斷垂直向上射出粒子，噴射口面積為10 m×10 m。假設風向沿x軸正方向，且y方向速度梯度為0，100 m以上風速為參考風速，100 m以下風速隨高度z的變化滿足大氣邊界層風廓線形式，即

(12)

式(12)中：u表示z高度處的風速，m/s；u*表示摩擦速度，m/s，u*的取值可以根據100 m高度處參考風速求出；K為馮卡門常數，取值為0.41；z0=0.1，代表粗糙度長度，m。

圖1 計算區域設置Fig.1 Computational domain configurations

本次模擬共設置a、b、c三組算例進行對比，其入口參考風速分別為2、5和10 m/s，以討論不同平均風速對擴散的影響。模式模擬的基本參數設置如表1所示。整個計算域網格數共計500×40×100=200 000個，模擬時長為1 200 s，噴射口粒子的噴射速度為垂直向上2 m/s，每秒噴射粒子1 000個，粒徑大小為0～200 μm，平均粒徑為100 μm，且假設其滿足標準正態分布。地面反射系數設為0.8，入口施加隨機脈動(±0.02, ±0.02,±0.01)m/s。計算中所選取的大渦模型為Smagorinsky模型。為保證數值穩定及時間計算精度，最大庫朗數設為0.6，最大時間步為1 s，并且每10 s輸出一次計算結果。

2 計算結果與分析

選取流場穩定狀態1 200 s時的數據進行分析，輸出11 m高度風速圖以體現流線變化情況。選取11 m高度處流線的原因是該處流場可充分體現源項上方的氣流變化情況。圖2所示為11 m高度處

表1 模式基本輸入參數

圖2 流場穩定時(t=1 200 s)流線分布Fig.2 The distribution of the stream line when the flow field becomes steady (t=1 200 s)

圖3 算例a中t=1 200 s時下風向不同距離截面處污染物的數濃度空間分布Fig.3 The spatial distribution of pollutant number concentration at t=1 200 s in case a

流線分布示意圖，圖中可以明顯看到：由于污染源對入口風速存在阻擋作用，源項周圍出現明顯繞流作用，源項后方出現旋轉方向相反的對稱渦旋。該繞流現象會導致污染物擴散時的分裂現象，具體變化情況在后文中詳細描述。

圖3所示為算例a(入口風速為2 m/s)t=1 200 s時刻計算域下風向1 000、1 100、1 200、1 300 m處污染物的數濃度空間分布。由于y方向速度梯度為零，且入口風速較小，污染物始終分裂成兩股，未發生匯合，兩股污染物之間距離為5～20 m。距離入口1 000 m處污染物數濃度最大值為78個/m3，1 100 m處數濃度最大值為70個/m3，1 200 m處數濃度最大值為80個/m3，1 300 m處污染物數濃度最大值為87個/m3。

圖4所示為算例b(風速大小為5 m/s)t=1 200 s時刻計算域下風向不同距離處污染物濃度空間分布。圖中可見：隨著風速的增大，分裂的兩股污染物之間的距離減小，主體污染物距離約在2～10 m左右。下風向1 000 m處污染物數濃度最大值為300個/m3，1 100 m處污染物數濃度最大值為212個/m3，1 200 m處數濃度最大值為227個/m3，1 300 m處污染物數濃度最大值為257個/m3。

圖5所示為算例c(風速大小為10 m/s)計算時間達到1 200 s時計算域內下風向不同距離污染物數濃度空間分布。從圖5中可以看到：風速的增大使得噴射口對于入口風的阻擋作用減小，該算例中污染物在擴散過程中未發生分裂；算例c中下風向1 000 m處污染物數濃度最大值為413個/m3, 1 100 m處污染物數濃度最大值為502個/m3, 1 200 m處污染物數濃度最大值為407個/m3, 1 300 m處污染物數濃度最大值為409個/m3。這些算例中污染物濃度隨風速增大而增大的原因在于入口風速的增大使得污染物粒子在計算域內停留時間變短，從而導致粒子位移距離減小，濃度增大。

圖5 算例c中t=1 200 s時下風向不同距離截面處污染物的數濃度空間分布Fig.5 The spatial distribution of pollutant number concentration at t=1 200 s in case c

圖6 污染物中心高度-時間變化曲線Fig.6 The temporal change of the height of the pollutant center

圖6顯示了下風向1 500 m處污染物數濃度最大值的高度隨時間的變化趨勢。圖中可見，隨著湍流的擾動作用，算例a污染物中心(濃度最大值)高度在63～76 m，算例b污染物中心高度約在46～52 m，算例c污染物中心高度在36～41 m。因此，風速的增大使得污染物濃度最大值高度減小，且波動范圍減小。圖7所示為下風向1 500 m處不同風速條件下粒子最大抬升高度隨時間變化曲線。隨著風速的增大，粒子最大抬升高度逐漸減小。2 m/s風速的情況下，粒子最大抬升高度可達95 m，5 m/s風速下粒子最大抬升高度為59 m，10 m/s風速下粒子最大抬升高度為41 m，且風速較大時最大抬升高度變化趨勢較為穩定。將污染物中心濃度值衰減至5%時的距離定義為影響半徑，圖8顯示了1 500 m處污染物最大影響半徑隨時間的變化情況。從圖中可以看到，2 m/s時影響半徑較大，數值在17～27 m，最大可達27 m；而5 m/s風速時影響半徑最大為12 m，10 m/s風速時最大影響半徑為7 m。風速的增大，使得影響半徑逐漸減小。

圖7 顆粒最大抬升高度-時間變化曲線 Fig.7 The temporal change of the maximum height that can be reached by the elevated particles

圖8 污染物影響半徑隨時間變化曲線Fig.8 The temporal change of the influencing radius of the pollutants

綜合圖6、圖7、圖8分析發現，在入口湍流條件不變的情況下，風速的增大會減小污染物粒子在計算域內的停留時間，由于y方向未設置速度梯度，污染物縱向擴散能力減弱，污染源中心濃度值增大，污染物抬升高度降低，污染物擴散范圍較小，即風速增大時，污染物軸向(最大影響半徑)徑向(抬升高度)擴散都呈現減小的趨勢。這一模擬結果，與文獻[14]中對單點源煙囪排放的顆粒物擴散的模擬結果相一致。

3 結論

應用MP-PIC方法模擬理想條件下不同風速條件對單點源污染源的污染物濃度時空變化、粒子最大抬升高度、污染物擴散范圍的影響，得到如下結論。

(1)MP-PIC方法可以極大地提高計算效率，可以應用于較大尺度的大氣污染物擴散模擬。

(2)CFD模擬中污染源對入口風的阻擋作用十分重要。低風速下，受污染源阻擋作用，污染源周邊形成方柱繞流現象，污染物在下風方向分裂成為兩股，污染物因湍流作用在計算域內波動。2 m/s風速時，分裂的兩股污染物之間距離在5～20 m，5 m/s風速時，分裂的兩股污染物距離為2～10 m，而風速達到10 m/s分裂狀況消失。

(3)風速的增大會減少污染物在計算域內的停留時間，使得粒子移動距離減小，從而導致污染源中心濃度值增大，污染物抬升高度與擴散范圍降低。即風速增大時，污染物軸向(最大影響半徑)徑向(抬升高度)擴散都呈現減小的趨勢。