徐路程，肖凱濤，宋偉偉

?

CFD方法在刺激劑擴(kuò)散模擬中的有效性驗證

徐路程，肖凱濤，宋偉偉

（防化研究院，北京，102205）

基于計算流體力學(xué)（CFD）方法，對某型刺激劑在平坦開闊地區(qū)的擴(kuò)散進(jìn)行模擬，并與實驗結(jié)果進(jìn)行對比，從而驗證了CFD方法對于模擬刺激劑在各種風(fēng)速條件下擴(kuò)散的有效性，為其在復(fù)雜地形環(huán)境中的模擬奠定理論基礎(chǔ)。

刺激劑；擴(kuò)散模擬；計算流體力學(xué)；驗證

刺激劑可以通過燃燒、爆炸、機(jī)械噴灑等方式施放到大氣中，在大氣中擴(kuò)散形成氣溶膠，由于在復(fù)雜環(huán)境中難以對刺激劑進(jìn)行測試試驗，因此很難了解不同的氣象條件下、復(fù)雜環(huán)境中刺激劑能夠達(dá)到的效能。而數(shù)值模擬方法能夠模擬出刺激劑在不同地形、不同氣象條件下的擴(kuò)散，進(jìn)而為其在復(fù)雜環(huán)境中的使用提供技術(shù)支持。

目前，國內(nèi)外已有多項研究使用計算流體力學(xué)方法對氣溶膠污染物的擴(kuò)散進(jìn)行模擬[1-4]，這些模擬都是在穩(wěn)態(tài)條件（即假定各物理量隨時間基本保持穩(wěn)定）下進(jìn)行的，顯然不適合模擬刺激劑的施放過程的瞬態(tài)性。因此，本文以某型刺激劑為例，使用計算流體力學(xué)方法對平坦開闊地條件下的刺激劑的施放進(jìn)行模擬，與實驗結(jié)果進(jìn)行對比，以驗證計算流體力學(xué)方法對于刺激劑擴(kuò)散模擬的適用性，并為這種方法運用于復(fù)雜地形環(huán)境下刺激劑擴(kuò)散模擬奠定理論基礎(chǔ)。

1 計算流體力學(xué)方法

1.1 幾何模型

用于數(shù)值模擬的幾何模型如圖1所示，計算域為140m×140m的正方形區(qū)域，在圖中釋放點處建立直角坐標(biāo)系，距離計算區(qū)域的右邊界和下邊界分別為10m，模擬的風(fēng)向為東偏南25°。由于氣溶膠擴(kuò)散在空間中是一個三維的過程，并且擴(kuò)散過程一般都發(fā)生在大氣邊界層以內(nèi)。因此，本模型中考慮垂直方向的計算高度為40m，于是，三維的計算區(qū)域為140m×140m×40m。

圖1 計算區(qū)域示意圖

1.2 數(shù)學(xué)模型

數(shù)學(xué)模型基于以下假定：（1）流動為不可壓縮流動。由于大氣邊界層中的流速遠(yuǎn)小于聲速，因此可以假定流動為不可壓縮流動；（2）流動為穩(wěn)態(tài)流動。文獻(xiàn)[5]中采用了這種假定，這種穩(wěn)態(tài)實際上是一種“偽穩(wěn)態(tài)”（pseudo- steady-state），能夠代表邊界層中一段時間內(nèi)相對穩(wěn)定的一種流動狀態(tài)，這種狀態(tài)也是適合于刺激劑擴(kuò)散的一種風(fēng)場條件。

湍流是大氣邊界層中的主要運動形態(tài)，對動量輸運、熱量輸送以及物質(zhì)輸送起著主要作用[6]。而計算流體力學(xué)方法中目前有3類方法模擬湍流[7]：直接數(shù)值模擬（DNS，Direct Numerical Simulation）方法，大渦模擬（LES，Large Eddy Simulation）方法，雷諾時均方程模擬（RANS，Reynolds Average Numerical Simulation）方法。這3種方法的區(qū)別在于模擬湍渦尺度范圍的不同。RANS方法從實際工程應(yīng)用的角度出發(fā)，將物理量都分解為平均項和脈動項，但由于引入了脈動項，就需要引入相應(yīng)的方程來描述脈動項，也就是要考慮加入所謂的湍流模型來使方程組封閉。由于較高的穩(wěn)定性和效率，-模型及其變種是最為常用的湍流模型[8]。本研究中將采用RANS方法中的-模型對湍流進(jìn)行模擬。

基于以上假定和湍流模型的選取，流體的控制方程組（張量形式）如下：

（2）

1.3 邊界條件

邊界條件的設(shè)定如圖2所示。

（1）東側(cè)、南側(cè)邊界設(shè)定為速度邊界條件，并且采用《環(huán)境影響評價技術(shù)導(dǎo)則大氣環(huán)境部分》所建議的冪指數(shù)風(fēng)廓線，在計算中通過使用Fluent軟件的UDF（User Defined Function，用戶自定義函數(shù)）實現(xiàn)：

實驗中風(fēng)速一般是由距地2m高的風(fēng)速儀得到的，而上述公式中需要的是10m高度的風(fēng)速。因此，需要首先推導(dǎo)2m風(fēng)速和10m風(fēng)速的關(guān)系：

10≈1.252 72（6）

模擬中2m高度的風(fēng)速的取值以及對應(yīng)的10m高處的風(fēng)速大小如表1所示。

圖2 三維幾何模型及邊界條件

表1 2m高度風(fēng)速和10m高度風(fēng)速的對應(yīng)關(guān)系 (m·s-1)

Tab.1 Correspondence table between2and10

模擬的風(fēng)向取為實驗時較為穩(wěn)定的風(fēng)向：東偏南25°，即與軸負(fù)向成25°，如圖1所示。溫度梯度取為-0.3K/m，2m高度處測得的溫度為8℃，因此垂直方向的溫度分布為：

()=281.15+0.3（-2） （7）

溫度廓線在Fluent中同樣使用UDF實現(xiàn)。

（2）西側(cè)、北側(cè)邊界設(shè)定為外流邊界，即所有的物理量在此邊界條件處都滿足梯度為零的條件。

（3）上邊界設(shè)定為對稱邊界，即此處物理量的法向速度為零，法向梯度也為零。

（4）下邊界設(shè)定為壁面邊界，即無滑移邊界，這種邊界代表了在邊界上的法向和切向的速度均為零。而根據(jù)溫度分布的要求，此處的溫度應(yīng)設(shè)置為：280.55K。

1.4 網(wǎng)格劃分

本研究使用Gambit軟件對幾何模型進(jìn)行網(wǎng)格劃分，Gambit提供了多種網(wǎng)格劃分的方式，而由于要在釋放點附近進(jìn)行局部的加密，因此選擇能夠靈活劃分的非結(jié)構(gòu)網(wǎng)格。本研究中采用如下網(wǎng)格劃分形式：水平方向上，靠近釋放點處的網(wǎng)格大小設(shè)置為0.5m，遠(yuǎn)離釋放點處的網(wǎng)格大小設(shè)置為2m，這是因為釋放點周圍的物理量變化會較快，較密的網(wǎng)格能夠較為準(zhǔn)確地表現(xiàn)出這種變化；垂直方向上由于在靠近壁面處的物理量變化較快，因此要布置稍密的網(wǎng)格，而在遠(yuǎn)離下邊界的上層，可以布置較為稀疏的網(wǎng)格。于是，在靠近地面處布置初始網(wǎng)格的尺寸為0.5m，垂直方向上共布置40個網(wǎng)格，沿著垂直向上的方向均勻增長；整個三維空間區(qū)域使用非規(guī)則的三棱柱網(wǎng)格進(jìn)行劃分，整個計算區(qū)域共108×104個網(wǎng)格。

1.5 流場求解方法

本研究使用Ansys Fluent 12.0對控制方程進(jìn)行求解。（1）控制方程離散：動量方程、能量方程、湍動能方程和湍流耗散率方程使用二階迎風(fēng)格式進(jìn)行離散，具有較好的精確度；壓力使用body-force - weighted方法，這種離散方法適用于存在自然對流的情況。（2）壓力速度耦合方式：選擇SIMPLEC算法對壓力速度進(jìn)行耦合[7]。（3）殘差設(shè)定：為了保證計算結(jié)果的收斂性，將能量方程的收斂準(zhǔn)則設(shè)定為10-6（標(biāo)準(zhǔn)化殘差），其他方程的收斂準(zhǔn)則設(shè)定為10-3（標(biāo)準(zhǔn)化殘差），同時，在計算過程中還對空間中一點進(jìn)行監(jiān)測，2m高度風(fēng)速大小為1m/s時的監(jiān)測結(jié)果如圖3所示，其他風(fēng)速條件下的風(fēng)場收斂監(jiān)測結(jié)果是類似的。從標(biāo)準(zhǔn)化殘差圖中可以看到各個指標(biāo)基本收斂，但收斂到什么程度并不是很直觀。收斂的最終目的是使計算域中的物理量無限地逼近真解，因此，采用對一點監(jiān)測物理量的辦法能夠更加清楚的說明這一點，如圖4所示，計算域中空間坐標(biāo)為(0,0,2)的點的速度大小已經(jīng)在迭代的過程中逐漸收斂。

圖3 標(biāo)準(zhǔn)化殘差收斂圖（u2=1m/s）

圖4 流場中一點速度大小監(jiān)測圖（u2=1m/s）

圖5 跨風(fēng)截面位置示意圖

圖6 跨風(fēng)截面風(fēng)速矢量圖（u2=1m/s）

圖5為跨風(fēng)截面位置示意圖，從跨風(fēng)截面風(fēng)速矢量圖6可以看到，在整個跨風(fēng)截面上，風(fēng)速分布基本能夠保持與入口相同的穩(wěn)定狀態(tài)，即整個計算區(qū)域的風(fēng)廓線基本上保持與入口相同。這里與上相同，只使用2=1m/s的情況為例進(jìn)行說明。

1.6 氣溶膠擴(kuò)散模擬

拉格朗日方法相對精度較高，求解復(fù)雜程度較低，而且能夠適用于復(fù)雜地形的情況。因此，本研究中選用這種方法來模擬刺激劑的擴(kuò)散。刺激劑的施放源的尺寸相比整個計算域的尺寸可以忽略不計，因此研究中考慮使用點源來模擬釋放源，點源的屬性如下：施放位置位于空間原點上方0.1m高度位置，空間坐標(biāo)為(0,0,0.1)；施放方向為沿著風(fēng)速方向，空間向量為(-0.906,0.423,0)；施放時間為0~30s；粒子粒徑為10μm；燃燒粒子的溫度假定為100℃（373.15K）；初始速度相比風(fēng)速較小，假定為0.1m/s；錐形點源的半角為15°，出口半徑為0.01m；防暴彈總藥量為50g，燃燒后有效作用率為21％，即有效藥量為13.5g，由此可以估算得到：彈藥的平均質(zhì)量流率為0.000 45kg/s；彈藥的材料密度為1 296kg/m3。模擬擴(kuò)散的總時間為60s，時間步長取為0.1s，共600個時間步?？梢允褂每臻g中一點的1min劑量（mg·min/m3）[10]來衡量刺激劑在1min內(nèi)的累計作用效果：

式（8）中的約等號是表示用數(shù)值積分來近似表示解析積分，具體的數(shù)學(xué)含義是：通過計算可以得到0.1s到60s，600個時刻的濃度計算結(jié)果，可以記為c(=1,…,600)，可以假定用第個時刻的瞬時濃度代表從-1時刻到時刻之間的平均濃度，于是1min內(nèi)的平均濃度可以表示為（8）最右端所示。

而對于關(guān)注的=1.5m高的平面上的每個結(jié)點都可以通過式（8）求出1min劑量，但Fluent軟件沒有提供這一功能。因此，本研究中使用C#語言編寫程序，分別讀取各個時刻的計算結(jié)果文件，將平面上每個點的濃度數(shù)據(jù)累加，最后得到1min劑量。

2 實驗驗證

2.1 實驗材料

驗證實驗中所需要的測量儀器及測試彈藥的數(shù)量及用途如表2所示。

表2 實驗儀器及彈藥

Tab.2 Experimental apparatus and ammunition

2.2 實驗方法

實驗中，采樣器布設(shè)方案如圖7所示。實驗時，當(dāng)?shù)貧庀髼l件為：風(fēng)向東偏南25°，溫度8℃，溫度梯度約為0.3K/m。風(fēng)速與模擬中的2m高度的風(fēng)速一致。當(dāng)通過綜合氣象觀測系統(tǒng)測量到的風(fēng)速基本穩(wěn)定于期望風(fēng)速時，引爆測試彈，同時采樣器開始同步采樣。彈爆1min后，采樣器停止采樣。

圖7 標(biāo)場及采樣點布置方案示意圖

2.3 結(jié)果分析

本研究中涉及的刺激劑的有效1min劑量為0.86mg·min/m3。在實際使用中，決策者更為關(guān)注的是有效濃度能夠達(dá)到的空間范圍。因此，對1.5m高度平面（人體站立呼吸平面）上有效濃度以上的區(qū)域的長、寬、面積及最大劑量值進(jìn)行統(tǒng)計。

首先繪制閾值劑量的等值線以確定有效區(qū)域，閾值等值線的繪制由軟件Surfer 8.0實現(xiàn)。由于Fluent計算得到的結(jié)果不是規(guī)則的網(wǎng)格分布，而Surfer中使用的等值線生成算法是基于結(jié)構(gòu)化網(wǎng)格的。因此，在Surfer軟件中首先要將非結(jié)構(gòu)化的數(shù)據(jù)插值到結(jié)構(gòu)網(wǎng)格上，將非結(jié)構(gòu)數(shù)據(jù)格式轉(zhuǎn)換為結(jié)構(gòu)化的網(wǎng)格文件形式（.grd）?？死锝鸱ㄊ菍臻g數(shù)據(jù)求最優(yōu)、線性、無偏內(nèi)插估計量的方法，所繪制的濃度等值線既能保持?jǐn)?shù)據(jù)的細(xì)節(jié)又可盡力地表現(xiàn)數(shù)據(jù)的變化趨勢，能夠最大程度地利用所給的信息分析數(shù)據(jù)空間的邊界，具有很高的空間外推精度。因此，本研究中選擇克里金法對Fluent計算的數(shù)據(jù)進(jìn)行網(wǎng)格化處理。實驗數(shù)據(jù)的處理方法是類似的，不過，在網(wǎng)格化數(shù)據(jù)之前，要首先把極坐標(biāo)系表示的采樣點坐標(biāo)轉(zhuǎn)化為直角坐標(biāo)系表示的坐標(biāo)。Surfer對Fluent模擬結(jié)果和實驗數(shù)據(jù)的統(tǒng)計結(jié)果的對比曲線如圖8所示。

圖8 有效區(qū)域面積與風(fēng)速關(guān)系圖

對比圖8幾條曲線可以看到，各種指標(biāo)在各種風(fēng)速條件下都比較接近：有效區(qū)域長度的最大相對誤差為21.850％；有效區(qū)域?qū)挾鹊淖畲笙鄬φ`差為5.85％；有效區(qū)域面積的最大相對誤差為30.398％；有效區(qū)域的最大劑量的最大相對誤差為12.413％。

另外，數(shù)值模擬結(jié)果和實驗結(jié)果在隨風(fēng)速變化的規(guī)律上都體現(xiàn)出相同的趨勢。由圖8（a）可見：當(dāng)風(fēng)速小于4m/s時，隨風(fēng)速增加有效區(qū)域長度增加；風(fēng)速超過4m/s時，有效區(qū)域長度隨風(fēng)速增加減小。由圖8（b）可見：當(dāng)風(fēng)速小于2.5m/s時，隨風(fēng)速增加，有效區(qū)域?qū)挾仍黾樱伙L(fēng)速超過2.5m/s時，有效區(qū)域?qū)挾入S風(fēng)速增加減小。由圖8（c）可見：當(dāng)風(fēng)速小于3.5m/s時，隨風(fēng)速增高，有效區(qū)域面積變大；風(fēng)速超過3.5m/s時，有效區(qū)域面積隨風(fēng)速增高而變小。由圖8（d）可見：對于所有風(fēng)速，隨著風(fēng)速增加，最大劑量值會一致性地下降。

3 結(jié)語

本文首次將計算流體力學(xué)方法應(yīng)用于刺激劑的擴(kuò)散模擬，并將模擬結(jié)果與實驗結(jié)果進(jìn)行了對比。對比顯示在各個風(fēng)速條件下，模擬結(jié)果與實驗結(jié)果的各指標(biāo)的誤差都不超過40％，這證明了計算流體力學(xué)方法用于刺激劑擴(kuò)散模擬的可行性和有效性，從而表明計算流體力學(xué)方法可以作為復(fù)雜地形條件下刺激劑作用效能評估研究的理論基礎(chǔ)。

[1] Riddle A, Carruthers D, Sharpe A, et al. Comparisons between FLUENT and ADMS for atmospheric dispersion modelling[J]. Atmospheric Environment, 2004, 38(7): 1 029-1 038.

[2] Di Sabatino S, Buccolieri R, Pulvirenti B, et al. Flow and pollutant dispersion in street canyons using FLUENT and ADMS-Urban[J]. Environmental Modeling & Assessment, 2008, 13(3): 369-381.

[3] Pullen J, Boris J P, Young T, et al. A comparison of contaminant plume statistics from a Gaussian puff and urban CFD model for two large cities[J]. Atmospheric Environment, 2005, 39(6): 1 049-1 068.

[4] 金穎, 周偉國. 煙氣擴(kuò)散的 CFD 數(shù)值模擬[J]. 安全與環(huán)境學(xué)報, 2002, 2(1): 21-23.

[5] Cheng W C, Liu C H, Leung D Y C. On the correlation of air and pollutant exchange for street canyons in combined wind-buoyancy-driven flow[J]. Atmospheric Environment, 2009, 43(24):3 682-3 690.

[6] 盛裴軒,毛節(jié)泰，李建國,等.大氣物理學(xué)[M].北京：北京大學(xué)出版社,2003.

[7] 吳清松.計算熱物理引論[M].合肥:中國科學(xué)技術(shù)大學(xué)出版社,2009.

[8] Li X X, Liu C H, Leung D Y C, et al. Recent progress in CFD modelling of wind field and pollutant transport in street canyons[J]. Atmospheric Environment, 2006, 40(29): 5 640- 5 658.

[9] 蔣維楣,孫鑒寧,曹文俊.空氣污染氣象學(xué)教程（第二版）[M].北京:氣象出版社,2004.

[10] 《核生化防護(hù)大辭典》編輯部.核化生防護(hù)大辭典[M].上海：上海辭書出版社,2000.

The Validation of CFD Method in Simulation of Irritant Dispersion

XU Lu-cheng, XIAO Kai-tao, SONG Wei-wei

(Research Institution of Chemical Defense, Beijing，102205)

Based on the principle of computational fluid dynamics (CFD) method, the dispersion of some kind of irritant in flat, open terrain was simulated. By comparing with experimental result, the effectiveness of CFD method in modeling irritant dispersion under the condition of different wind velocities is verified, which establishes the theoretical foundation for its simulation in complex terrain.

Irritant；Dispersion simulation；Computational fluid dynamics；Validation

1003-1480（2014）05-0042-05

TQ567.5

A

2014-05-13

徐路程（1990-），男，在讀碩士研究生，主要從事防暴彈藥性能評價及使用技術(shù)研究。

總裝備“十二五”部預(yù)先研究課題。