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

徐路程，肖凱濤，宋偉偉

?

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

徐路程，肖凱濤，宋偉偉

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

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

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

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

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

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

1.1 幾何模型

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

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

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

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

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

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

（2）

1.3 邊界條件

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

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

實(shí)驗(yà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)速的取值以及對(duì)應(yīng)的10m高處的風(fēng)速大小如表1所示。

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

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

Tab.1 Correspondence table between2and10

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

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

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

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

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

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

1.4 網(wǎng)格劃分

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

1.5 流場(chǎng)求解方法

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

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

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

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

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

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

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

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

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

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

2 實(shí)驗(yàn)驗(yàn)證

2.1 實(shí)驗(yàn)材料

驗(yàn)證實(shí)驗(yàn)中所需要的測(cè)量?jī)x器及測(cè)試彈藥的數(shù)量及用途如表2所示。

表2 實(shí)驗(yàn)儀器及彈藥

Tab.2 Experimental apparatus and ammunition

2.2 實(shí)驗(yàn)方法

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

圖7 標(biāo)場(chǎng)及采樣點(diǎn)布置方案示意圖

2.3 結(jié)果分析

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

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

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

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

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

3 結(jié)語(yǔ)

本文首次將計(jì)算流體力學(xué)方法應(yīng)用于刺激劑的擴(kuò)散模擬，并將模擬結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比。對(duì)比顯示在各個(gè)風(fēng)速條件下，模擬結(jié)果與實(shí)驗(yàn)結(jié)果的各指標(biāo)的誤差都不超過(guò)40％，這證明了計(jì)算流體力學(xué)方法用于刺激劑擴(kuò)散模擬的可行性和有效性，從而表明計(jì)算流體力學(xué)方法可以作為復(fù)雜地形條件下刺激劑作用效能評(píng)估研究的理論基礎(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] 金穎, 周偉國(guó). 煙氣擴(kuò)散的 CFD 數(shù)值模擬[J]. 安全與環(huán)境學(xué)報(bào), 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é)泰，李建國(guó),等.大氣物理學(xué)[M].北京：北京大學(xué)出版社,2003.

[7] 吳清松.計(jì)算熱物理引論[M].合肥:中國(guó)科學(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-），男，在讀碩士研究生，主要從事防暴彈藥性能評(píng)價(jià)及使用技術(shù)研究。

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