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

0 引言

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

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

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

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

本文考慮到在實(shí)際情況中污染源實(shí)時(shí)定位對(duì)時(shí)間的要求，研究粗網(wǎng)格在污染源定位中的準(zhǔn)確性。網(wǎng)格數(shù)量的減少會(huì)引起流場(chǎng)計(jì)算的誤差增大，但聯(lián)合概率密度法對(duì)氣流場(chǎng)數(shù)據(jù)的精度要求不是很高，因此本文主要驗(yàn)證了流場(chǎng)精確度下降后的定位準(zhǔn)確性，從而在實(shí)際工程中可以減少污染物尋源的時(shí)間，實(shí)現(xiàn)實(shí)時(shí)定位。

1 研究原理

1.1 控制方程介紹

本研究中污染源定位模型主要是通過(guò)CFD粗網(wǎng)格模型數(shù)值模擬的方法驗(yàn)證其有效性。

控制方程的通用形式：

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

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

div(ρuφ)——對(duì)流項(xiàng)；

di (Γgradφ)——擴(kuò)散項(xiàng)；

?！獜V義擴(kuò)散系數(shù)；

S——廣義源項(xiàng)。

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

1.2 污染物尋源反計(jì)算

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

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

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

τ——逆向時(shí)間；

Vj,τ——表示某一時(shí)刻τ的氣流場(chǎng)；

vc——污染源在空氣中的有效擴(kuò)散系數(shù)；

δ（x）——脈沖函數(shù)；

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

式中 N——探頭處監(jiān)測(cè)到的濃度數(shù)據(jù)的個(gè)數(shù)；

τ0——污染源的釋放時(shí)間；

M0——瞬時(shí)污染源所釋放的污染物質(zhì)量總量；

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

通過(guò)求解式（4）、（5）可計(jì)算得污染源的概率聯(lián)合分布（CALP）。

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

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

圖1 基于CFD污染源定位反計(jì)算流程

2 案例驗(yàn)證

本文選用一個(gè)簡(jiǎn)化的二維辦公室模型來(lái)驗(yàn)證上述的理論。由于辦公室空間相對(duì)封閉，假設(shè)有污染源存在，若不能快速定位，后果將不堪設(shè)想。在此前的研究中，就已經(jīng)在網(wǎng)格無(wú)關(guān)解的模型下使用概率反計(jì)算法精確的定位辦公室污染源的位置，但耗費(fèi)大量的時(shí)間，在實(shí)際中對(duì)于人員的安全十分不利。因此，在辦公室模型網(wǎng)格無(wú)關(guān)解的基礎(chǔ)上減少網(wǎng)格的數(shù)量，建立快速定位的粗網(wǎng)格模型。

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

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

圖2 辦公室模型在測(cè)線A處不同網(wǎng)格下的CFD數(shù)值模擬速度結(jié)果

2.1 網(wǎng)格數(shù)量的選取

針對(duì)該模型，本文分別設(shè)置了5種不同數(shù)量的網(wǎng)格，分別為 15×15，30×30，50×50，80×80 以及100×100的網(wǎng)格，首先進(jìn)行網(wǎng)格無(wú)關(guān)解的計(jì)算，然后以此為基礎(chǔ)確定粗網(wǎng)格模型。在辦公室模型的中間位置（5 m處）設(shè)置了一條測(cè)線A，在測(cè)線A處記錄不同時(shí)刻氣流的大小。這5種網(wǎng)格在測(cè)線A處的計(jì)算結(jié)果如圖2所示。

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

2.2 粗網(wǎng)格模型污染物尋源的準(zhǔn)確性

通過(guò)上述辦公室模型中測(cè)線A處的不同時(shí)刻氣流場(chǎng)大小的分析，發(fā)現(xiàn)30×30的網(wǎng)格相對(duì)于網(wǎng)格無(wú)關(guān)解，可以粗略地描述其速度場(chǎng)。對(duì)不同網(wǎng)格數(shù)獲得的氣流場(chǎng)進(jìn)行尋源反計(jì)算，并比較30×30的網(wǎng)格和網(wǎng)格無(wú)關(guān)解的污染源定位結(jié)果，驗(yàn)證粗網(wǎng)格模型對(duì)污染源定位的準(zhǔn)確性。

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

圖3 探頭1多個(gè)污染物濃度數(shù)據(jù)聯(lián)合探頭2多個(gè)污染物濃度數(shù)據(jù)CALP定位結(jié)果

3 粗網(wǎng)格模型討論

通過(guò)以上粗網(wǎng)格的污染物尋源反計(jì)算，可以發(fā)現(xiàn)粗網(wǎng)格可以對(duì)污染源進(jìn)行較為精確的定位，對(duì)此，為了進(jìn)一步驗(yàn)證粗網(wǎng)格模型可以準(zhǔn)確定位污染源的位置，比較 15×15,30×30和80×80的網(wǎng)格模型的速度場(chǎng)。

已知80×80的網(wǎng)格為網(wǎng)格無(wú)關(guān)解，30×30的網(wǎng)格是所選取的粗網(wǎng)格，在網(wǎng)格無(wú)關(guān)解下通過(guò)CFD模擬得到的速度場(chǎng)能夠準(zhǔn)確地描述辦公室模型內(nèi)部流動(dòng)特征，以其為基準(zhǔn)，并以此來(lái)比較15×15網(wǎng)格和粗網(wǎng)格下的速度場(chǎng)，對(duì)比分析三者速度場(chǎng)主要渦流的大小和方向，從而驗(yàn)證粗網(wǎng)格是否可以準(zhǔn)確的定位真實(shí)污染源的位置。不同網(wǎng)格數(shù)的速度場(chǎng)如圖4～7所示。

圖4 t=60 s時(shí)刻速度云圖

圖5 t=120 s時(shí)刻速度云圖

圖6 t=180s 時(shí)刻速度云圖

圖7 t=240 s時(shí)刻速度云圖

通過(guò)圖4～7各個(gè)時(shí)刻的速度場(chǎng)云圖對(duì)比分析可知，15×15網(wǎng)格由于其尺寸相對(duì)較大，造成其速度場(chǎng)的主要渦流和網(wǎng)格無(wú)關(guān)解的差距比較明顯，尤其是在靠近進(jìn)、出風(fēng)口處；相比較而言，所選取的粗網(wǎng)格的主要渦流的大小和方向都和網(wǎng)格無(wú)關(guān)解的速度場(chǎng)的十分相似，沒(méi)有明顯差別，可以確定粗網(wǎng)格是可以對(duì)污染源進(jìn)行準(zhǔn)確的定位。

通過(guò)以上的分析，確定本文選取的粗網(wǎng)格可以準(zhǔn)確的對(duì)污染源進(jìn)行定位。最后，通過(guò)統(tǒng)計(jì)不同網(wǎng)格使用CFD模擬得到速度場(chǎng)所需要的運(yùn)行時(shí)間，結(jié)果見(jiàn)表1。

表1 不同網(wǎng)格下CFD模擬求解所需時(shí)間對(duì)比

計(jì)算采用的電腦型號(hào)為聯(lián)想ThinkPad SL410，處理器為英特爾酷睿2雙核T6570@2.1GHz筆記本處理器，內(nèi)存為4 GB（爾必達(dá)DDR3 1066 MHz）。

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

僅是正向速度場(chǎng)計(jì)算的時(shí)間，粗網(wǎng)格模型就相較于網(wǎng)格無(wú)關(guān)解節(jié)約了70%以上的時(shí)間，完整的反計(jì)算要在此基礎(chǔ)上進(jìn)行速度取反后，計(jì)算伴隨概率SALP和聯(lián)合概率CALP，要耗費(fèi)更長(zhǎng)的時(shí)間。

4 結(jié)語(yǔ)

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