基于粗網格CFD模擬方法的室內空氣污染物實時尋源反計算研究

（上海理工大學 環(huán)境與建筑學院，上海 200093）

0 引言

在現代社會中，人們在室內的時間達到了90%以上[1]，室內環(huán)境的空氣質量安全顯得尤為重要。污染物、火災以及危險化學品意外泄漏等事故發(fā)生時，濃煙和毒氣在短時間內擴散，此時需要快速定位污染源的位置，以保證人員的安全。

針對于如何快速、準確定位污染源位置的問題，反計算方法可以很好解決。通過分析探頭處的污染物濃度信息和相應的氣流組織，逆向推斷出污染源的位置和釋放時間?，F在反計算已被廣泛應用到相關實際問題中，比如將反計算應用于大氣污染源位置的辨識[2-5]等相關研究之中；火災發(fā)生時，通過相應的定位模型快速定位到火災源頭[6]。針對于反計算求解問題，Zhai將其歸結為三類方法：正向求解法、反向求解法和概率反計算法[7]。

正向法是一種試湊法，假設所有可能的污染源位置和散發(fā)情況，并針對每種情況采用正向模擬求解計算結果，與實際情況最為相似的結果所對應的源的情況即為污染源位置最可能的解[8-10]；該方法效率低，耗時長。反向法從最終污染物的濃度場開始，采用負時間步長對污染物的擴散進行反推，最終獲得污染源的位置[7]；該方法需要有最終濃度場所有的信息作為輸入條件才能求解，在實際中這些信息很難被監(jiān)測到，導致該方法使用的場合受限。概率反計算法則是采用概率分析估算某一事件作為概率目標函數，通過最大或最小化該目標函數來確定污染源的相關信息[11]，這種方法相比較于前兩種，在運算強度、輸入信息以及計算的穩(wěn)定性方面都有明顯的優(yōu)勢；Liu等已通過該方法成功定位到室內污染源的位置[12-13]。

然而，上述的方法雖然有著較好的精確性，但在實時性方面卻表現的不理想，主要是由于CFD模擬獲得氣流場的時間很長，特別是對于非穩(wěn)態(tài)的流動問題。人們開始研究如何提高模擬速度，一些學者提出使用簡單的湍流模型[14]、采用先進的模擬軟件和硬件以及使用粗網格模型[15]等方法。其中，網格密度和分布是影響CFD模擬時間的關鍵因素，對模擬的結果也有較大的影響[16]。在室內污染物尋源反計算的過程中，通常采用網格無關解進行模擬，以準確定位污染源的位置，而在實際情況中，受到計算機存儲量、計算速度的限制，網格數量不可能無限大[17]。有研究顯示，一個機艙模型在使用RANS湍流模型進行CFD計算的情況下，當網格數精密到可以準確模擬出機艙內非穩(wěn)態(tài)氣流組織形式的時候，使用k-ε模型所需時間大約為11 h，使用其他更為精密的模型則需要2倍甚至更久[18]。CFD模擬通過增大網格數量來對區(qū)域內的氣流場和污染物擴散進行較為精確的描述，耗費了大量的時間。在室內污染物尋源反計算中，計算時間過長對于迅速定位污染源位置十分不利，因此需要研究一種快速有效的方法。相關文獻中表明，網格劃分的數量越多計算精度就越高，但工作量就越大、仿真分析效率就越低，所以不能采用無限增大網格數的方法來提高計算精度[19]；在CFD模擬中粗網格計算得到的結果與網格無關解的相差不大，并且可以大大減少運算的時間，假設網格無關解為1個單位的計算時間，則粗網格的計算時間要小于1%[16]。因此應該在保證精度的情況下盡可能的減少網格數，以快速定位污染源的位置。

本文考慮到在實際情況中污染源實時定位對時間的要求，研究粗網格在污染源定位中的準確性。網格數量的減少會引起流場計算的誤差增大，但聯合概率密度法對氣流場數據的精度要求不是很高，因此本文主要驗證了流場精確度下降后的定位準確性，從而在實際工程中可以減少污染物尋源的時間，實現實時定位。

1 研究原理

1.1 控制方程介紹

本研究中污染源定位模型主要是通過CFD粗網格模型數值模擬的方法驗證其有效性。

控制方程的通用形式：

式中 ?(ρφ)/?t——瞬態(tài)項；

φ——通用變量，可以代表速度、溫度，或是其他待求解變量；

div(ρuφ)——對流項；

di (Γgradφ)——擴散項；

?！獜V義擴散系數；

S——廣義源項。

對于不同的控制方程，φ，Γ和S分別具有不同的形式。

1.2 污染物尋源反計算

污染物尋源反計算本質上是污染物傳播特性的研究，通過相應的程序將動態(tài)氣流場的方向取反，在探測器所處位置釋放瞬時污染物，污染物將通過方向取反后的動態(tài)氣流組織進行傳播，在相應的時間內傳播到真實污染源的位置，然后依據污染物正向傳播時監(jiān)測到的污染物濃度數據，結合概率分布理論求取污染源位置在整個空間區(qū)域的概率分布，其中概率值最大的點所對應的位置即認為是污染源位置[20]。

本文采用概率反計算法，以尋找源的位置展開研究。在矩形區(qū)域中的任何位置只要有污染物存在都可以認為其是潛在的污染源的位置[21]。概率反計算的表達式如下所示：

式中ψ＊—,—伴隨概率因子；

τ——逆向時間；

Vj,τ——表示某一時刻τ的氣流場；

vc——污染源在空氣中的有效擴散系數；

δ（x）——脈沖函數；

利用探頭監(jiān)測到的污染物歷史濃度信息來提高污染源定位的精確度；其具體的表達式如下：

式中 N——探頭處監(jiān)測到的濃度數據的個數；

τ0——污染源的釋放時間；

M0——瞬時污染源所釋放的污染物質量總量；

fx（X；τ0，Xi，τi）為第i個監(jiān)測探頭的SALP；為根據多個探測數據計算得到的基于污染物的釋放質量M0和位置處的概率分布。

通過求解式（4）、（5）可計算得污染源的概率聯合分布（CALP）。

為求解污染源位置的CALP，需要先求解得到污染源位置的SALP。求解SALP需要2個必要信息：非穩(wěn)態(tài)的速度場以及在室內環(huán)境中布置的探頭位置。其中，從式（3）中可以發(fā)現，式中的Vj,τ表示某一時刻的氣流場，本文中污染源所存在的非穩(wěn)態(tài)氣流場是隨時間變化的，因此網格的數量將影響計算的速度；采用網格無關解進行計算時，要耗費大量時間，若在保證尋源結果準確的前提下，在網格無關解的基礎上減少網格的數量得到粗網格，將會大大縮短尋源計算的時間。所需要的非穩(wěn)態(tài)的速度場、污染物擴散的濃度場均可通過CFD求解得到。

針對不同的具體案例，首先需要根據案例的具體參數進行建模，確定所設探頭的位置、數量及其他物理參數；隨后分別對網格無關解和粗網格模型進行CFD模擬，得到相應的動態(tài)氣流場；然后將動態(tài)氣流組織數據和探頭在不同時刻監(jiān)測得到的污染物濃度信息作為輸入信息，輸入反演計算程序。反演計算程序通過求解上文所闡述的公式，計算出污染源位置的伴隨概率（SALP）和聯合概率（CALP）。圖1示出基于CFD污染源定位反計算的流程。

圖1 基于CFD污染源定位反計算流程

2 案例驗證

本文選用一個簡化的二維辦公室模型來驗證上述的理論。由于辦公室空間相對封閉，假設有污染源存在，若不能快速定位，后果將不堪設想。在此前的研究中，就已經在網格無關解的模型下使用概率反計算法精確的定位辦公室污染源的位置，但耗費大量的時間，在實際中對于人員的安全十分不利。因此，在辦公室模型網格無關解的基礎上減少網格的數量，建立快速定位的粗網格模型。

辦公室的長為10 m，高為3 m，假設辦公室內有1人，1張桌子，桌子上放置1臺電腦，人坐在電腦前面辦公。其中假設人的發(fā)熱量為70 W，電腦的發(fā)熱量為200 W，在左墻的上方有一個水平條縫型送風口，送風速度0.1 m/s，送風溫度T=20 ℃。右墻的下方設有一個水平條縫型出風口。

在辦公室內安裝2個污染物監(jiān)測探頭：探頭1被安裝在天花板中間位置，探頭2被安裝在辦公室右上角。假設在窗戶下方有一個瞬時釋放的點污染源，源的大小為100單位。在t=0時刻開始釋放，與此同時空調開啟，送風口開始送風，在污染物的釋放過程中，空調的氣流組織還未達到穩(wěn)定狀態(tài)，時刻發(fā)生著變化，是一個典型的非穩(wěn)態(tài)氣流場。探測器從t=0時刻開始監(jiān)測其周圍環(huán)境污染物濃度的變化，直到整個過程結束。整個模擬過程持續(xù)240 s，分為8步，時間步長為30 s。采用Phoenics軟件，湍流模型選擇標準的k-ε模型。

圖2 辦公室模型在測線A處不同網格下的CFD數值模擬速度結果

2.1 網格數量的選取

針對該模型，本文分別設置了5種不同數量的網格，分別為 15×15，30×30，50×50，80×80 以及100×100的網格，首先進行網格無關解的計算，然后以此為基礎確定粗網格模型。在辦公室模型的中間位置（5 m處）設置了一條測線A，在測線A處記錄不同時刻氣流的大小。這5種網格在測線A處的計算結果如圖2所示。

從圖2比較5種網格數模擬結果，可以發(fā)現50×50網格、80×80網格和 100×100網格的差異較小，尤其是80×80網格和100×100網格的計算結果非常相似，可見80×80網格足以較好地描述該模型內部的流動特征，且在網格數增多的時候，模擬得到的非穩(wěn)態(tài)氣流場不再發(fā)生顯著的差異，因此可認為80×80的網格為網格無關解。

2.2 粗網格模型污染物尋源的準確性

通過上述辦公室模型中測線A處的不同時刻氣流場大小的分析，發(fā)現30×30的網格相對于網格無關解，可以粗略地描述其速度場。對不同網格數獲得的氣流場進行尋源反計算，并比較30×30的網格和網格無關解的污染源定位結果，驗證粗網格模型對污染源定位的準確性。

圖3示出探頭1多個污染物濃度數據聯合探頭2多個污染物濃度數據作為污染源定位的輸入部分，通過定位模型的反計算得到污染源的位置。通過對比圖3，發(fā)現兩者的定位結果差異很小，都可以十分準確地定位到真實污染源的位置。即使是在30×30的網格下，其預測可能是真實污染源的位置也集中在很小的區(qū)域，較為集中，準確性較高。在30×30的網格中，通過污染源定位模型計算得到的有可能是真實污染源的位置的最高概率也高達90%，非常接近真實污染源所處的位置，所以通過30×30的網格模型可以進行準確的定位。因此，選取30×30的網格為粗網格。

圖3 探頭1多個污染物濃度數據聯合探頭2多個污染物濃度數據CALP定位結果

3 粗網格模型討論

通過以上粗網格的污染物尋源反計算，可以發(fā)現粗網格可以對污染源進行較為精確的定位，對此，為了進一步驗證粗網格模型可以準確定位污染源的位置，比較 15×15,30×30和80×80的網格模型的速度場。

已知80×80的網格為網格無關解，30×30的網格是所選取的粗網格，在網格無關解下通過CFD模擬得到的速度場能夠準確地描述辦公室模型內部流動特征，以其為基準，并以此來比較15×15網格和粗網格下的速度場，對比分析三者速度場主要渦流的大小和方向，從而驗證粗網格是否可以準確的定位真實污染源的位置。不同網格數的速度場如圖4～7所示。

圖4 t=60 s時刻速度云圖

圖5 t=120 s時刻速度云圖

圖6 t=180s 時刻速度云圖

圖7 t=240 s時刻速度云圖

通過圖4～7各個時刻的速度場云圖對比分析可知，15×15網格由于其尺寸相對較大，造成其速度場的主要渦流和網格無關解的差距比較明顯，尤其是在靠近進、出風口處；相比較而言，所選取的粗網格的主要渦流的大小和方向都和網格無關解的速度場的十分相似，沒有明顯差別，可以確定粗網格是可以對污染源進行準確的定位。

通過以上的分析，確定本文選取的粗網格可以準確的對污染源進行定位。最后，通過統(tǒng)計不同網格使用CFD模擬得到速度場所需要的運行時間，結果見表1。

表1 不同網格下CFD模擬求解所需時間對比

計算采用的電腦型號為聯想ThinkPad SL410，處理器為英特爾酷睿2雙核T6570@2.1GHz筆記本處理器，內存為4 GB（爾必達DDR3 1066 MHz）。

從表1中可以發(fā)現，15×15網格下CFD模擬求解速度場所需要的時間最短，為5 min，但15×15的網格由于其網格尺寸較大，對細節(jié)的特征描述不夠，相較于網格無關解的速度場相差較大，不能準確定位到污染源的位置。80×80的網格無關解與30×30的粗網格污染源定位結果相似，都較為準確，雖然網格無關解描述的速度場更加接近真實情況，但是其求解速度場需要耗費大量的時間，不能實時求解；而對于30×30的粗網格，計算時間縮短的同時又可以十分準確地定位污染源的位置。

僅是正向速度場計算的時間，粗網格模型就相較于網格無關解節(jié)約了70%以上的時間，完整的反計算要在此基礎上進行速度取反后，計算伴隨概率SALP和聯合概率CALP，要耗費更長的時間。

4 結語

本文通過二維典型送風模型，利用不同網格數CFD模擬獲得的氣流場進行尋源反計算，以網格無關解的定位結果為基準研究粗網格尋源定位的準確性?？梢园l(fā)現，粗網格模型在非穩(wěn)態(tài)氣流場中可以準確的定位污染源的位置，而且僅與網格無關解的定位結果存在微小的誤差，可以忽略不計；使用粗網格模型進行CFD模擬相較于網格無關解節(jié)約了70%以上的時間，提高了計算速度，可在非穩(wěn)態(tài)氣流組織工況下快速定位到污染源的位置，為后續(xù)的救援贏得寶貴的時間。