徐路程，肖凱濤，宋偉偉

?

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

徐路程，肖凱濤，宋偉偉

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

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

刺激劑；擴散模擬；計算流體力學；驗證

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

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

1 計算流體力學方法

1.1 幾何模型

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

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

1.2 數(shù)學模型

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

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

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

（2）

1.3 邊界條件

邊界條件的設定如圖2所示。

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

實驗中風速一般是由距地2m高的風速儀得到的，而上述公式中需要的是10m高度的風速。因此，需要首先推導2m風速和10m風速的關系：

10≈1.252 72（6）

模擬中2m高度的風速的取值以及對應的10m高處的風速大小如表1所示。

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

表1 2m高度風速和10m高度風速的對應關系 (m·s-1)

Tab.1 Correspondence table between2and10

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

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

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

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

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

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

1.4 網(wǎng)格劃分

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

1.5 流場求解方法

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

圖3 標準化殘差收斂圖（u2=1m/s）

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

圖5 跨風截面位置示意圖

圖6 跨風截面風速矢量圖（u2=1m/s）

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

1.6 氣溶膠擴散模擬

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

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

而對于關注的=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 實驗方法

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

圖7 標場及采樣點布置方案示意圖

2.3 結(jié)果分析

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

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

圖8 有效區(qū)域面積與風速關系圖

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

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

3 結(jié)語

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

[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] 金穎, 周偉國. 煙氣擴散的 CFD 數(shù)值模擬[J]. 安全與環(huán)境學報, 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é)泰，李建國,等.大氣物理學[M].北京：北京大學出版社,2003.

[7] 吳清松.計算熱物理引論[M].合肥:中國科學技術大學出版社,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] 蔣維楣,孫鑒寧,曹文俊.空氣污染氣象學教程（第二版）[M].北京:氣象出版社,2004.

[10] 《核生化防護大辭典》編輯部.核化生防護大辭典[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-），男，在讀碩士研究生，主要從事防暴彈藥性能評價及使用技術研究。

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