基于OpenFOAM的化學(xué)危害擴(kuò)散預(yù)測求解器的開發(fā)與驗(yàn)證

韓朝帥， 諸雪征， 顧進(jìn)， 張宏遠(yuǎn)， 張赫， 左欽文

(1.陸軍防化學(xué)院， 北京 100000； 2.軍事科學(xué)院 防化研究院， 北京 100000；3.陸軍裝備部項(xiàng)目管理中心， 北京 100000)

0 引言

化學(xué)危害擴(kuò)散預(yù)測作為化學(xué)危害擴(kuò)散模擬的關(guān)鍵基礎(chǔ)，要將地形、障礙物、湍流等對危害擴(kuò)散的擾動(dòng)進(jìn)行精細(xì)描述，以求精準(zhǔn)可靠地預(yù)測未來時(shí)刻危害時(shí)空傳輸和擴(kuò)散情況，是國家核生化安全和國防安全領(lǐng)域的重要基礎(chǔ)性課題，對精確化提供防化保障輔助決策意義重大[1]。

化學(xué)危害擴(kuò)散具有擴(kuò)散尺度小、易受下墊面干擾，外部環(huán)境復(fù)雜、干擾因素多，數(shù)據(jù)來源雜、精度格式不一等特點(diǎn)，使得化學(xué)危害擴(kuò)散預(yù)測的復(fù)雜性和不確定性大大增強(qiáng)。Scargiali等[2]采用計(jì)算流體力學(xué)軟件CFX 4.4將重氣擴(kuò)散分解成污染物釋放前的流場模擬、污染物釋放時(shí)的擴(kuò)散模擬兩個(gè)步驟，開展了重氣擴(kuò)散精準(zhǔn)預(yù)測研究。Henry等[3]利用大氣模擬預(yù)測軟件AERMOD模型對微粒污染的輸運(yùn)和擴(kuò)散進(jìn)行模擬分析，驗(yàn)證了大氣擴(kuò)散模型(AERMOD)在模擬面源排放顆粒物輸運(yùn)方面的預(yù)測能力。Vianello等[4]采用計(jì)算流體力學(xué)軟件Fluent和工藝危險(xiǎn)源分析軟件(PHAST)開展城區(qū)氯氣排放的數(shù)值模擬研究，分析了兩種軟件的適用性。

然而，目前化學(xué)危害擴(kuò)散預(yù)測大多是通過商業(yè)軟件來進(jìn)行的，而現(xiàn)有商業(yè)軟件中基本沒有針對戰(zhàn)場化學(xué)危害擴(kuò)散模擬的功能，難以形成綜合性的態(tài)勢感知。在我國大力發(fā)展自主知識產(chǎn)權(quán)軟件的大背景下，計(jì)算流體力學(xué)軟件OpenFOAM的開源特性有助于打破軟件的版權(quán)壁壘，自主地開發(fā)研究工具。Deng等[5]將內(nèi)雙迭代高效算法代替interFOAM求解器中的PIMPLE算法，構(gòu)建適用于非穩(wěn)態(tài)兩相流問題求解的壓力速度耦合求解器，并通過算例驗(yàn)證了所建算法的正確性、收斂性和有效性。曾劍銳等[6]通過熱物性表制備與讀取、原熱物理派生類改寫、熱物性庫編譯3個(gè)步驟，構(gòu)建了適用于超臨界CO2流動(dòng)與傳熱問題的sCO2FOAM求解器，該研究對于OpenFOAM求解器的移植性和拓展性編譯的借鑒性較強(qiáng)。馬宇等[7]針對中子輸運(yùn)模擬的復(fù)雜性和耦合性，基于OpenFOAM軟件和有限體積法構(gòu)建了中子輸運(yùn)動(dòng)力學(xué)求解器ntkFOAM，實(shí)現(xiàn)了中子輸運(yùn)過程中傳熱傳質(zhì)的精細(xì)耦合計(jì)算。Abhishek等[8]將連續(xù)等離子體輸運(yùn)模型用于OpenFOAM軟件開發(fā)平臺(tái)，通過對幾何體并行模擬、網(wǎng)格剖分、動(dòng)量守恒模型進(jìn)行設(shè)計(jì)編譯，構(gòu)建SOMAFOAM求解器，實(shí)現(xiàn)了低溫等離子體在有限體積內(nèi)的連續(xù)擴(kuò)散模擬。

此外，化學(xué)危害擴(kuò)散預(yù)測還包括信息融合、源項(xiàng)反演、數(shù)據(jù)同化等3項(xiàng)關(guān)鍵技術(shù)，通過OpenFOAM軟件定制化開發(fā)化學(xué)危害擴(kuò)散預(yù)測求解器，對于研究這3項(xiàng)技術(shù)具有很強(qiáng)的互通拓展性。陳存楊等[9]在壓力- 隱式操作符分裂(PISO)FOAM求解器的基礎(chǔ)上引入濃度標(biāo)量場，開發(fā)用于重氣擴(kuò)散濃度預(yù)測的MyPISOFOAM求解器，但是其用于控制濃度標(biāo)量場的微分方程未考慮沉降、多源效應(yīng)，同時(shí)未植入相應(yīng)的信息融合、源項(xiàng)反演、數(shù)據(jù)同化等算法，導(dǎo)致預(yù)測精度不高。因此，本文在PISOFOAM瞬態(tài)求解器[10]的基礎(chǔ)上，將信息融合算法、四維變分同化算法(4DVAR)[11]、智能尋優(yōu)算法植入OpenFOAM軟件，詳細(xì)研究化學(xué)危害擴(kuò)散預(yù)測中基礎(chǔ)參數(shù)和數(shù)據(jù)導(dǎo)入并實(shí)現(xiàn)智能更新的方法，構(gòu)建用于化學(xué)危害擴(kuò)散預(yù)測的求解器ChdpFOAM，同時(shí)對新開發(fā)的求解器進(jìn)行應(yīng)用驗(yàn)證，以確保不會(huì)因參數(shù)設(shè)置不當(dāng)而導(dǎo)致預(yù)測失準(zhǔn)。

1 ChdpFOAM求解器

1.1 ChdpFOAM求解器工作流程設(shè)計(jì)

PISOFOAM瞬態(tài)求解器中，PISO算法通過迭代求解離散化的動(dòng)量預(yù)測方程和壓力泊松方程，實(shí)現(xiàn)速度與壓力的收斂耦合，進(jìn)而構(gòu)建擴(kuò)散空間流場分布[12]。在PISOFOAM瞬態(tài)求解器的基礎(chǔ)上，ChdpFOAM求解器綜合化學(xué)危害擴(kuò)散環(huán)境復(fù)雜、數(shù)據(jù)來源多維、邊界條件多變的特點(diǎn)，嵌入信息融合模塊、數(shù)據(jù)同化模塊、源項(xiàng)反演模塊，并通過智能尋優(yōu)算法實(shí)現(xiàn)邊界條件的智能更新和求解器的快速計(jì)算。ChdpFOAM求解器工作流程如圖1所示。

圖1 ChdpFOAM求解器工作流程Fig.1 Flow chart of ChdpFOAM solver

由圖1可知，在OpenFOAM軟件中構(gòu)建ChdpFOAM求解器的關(guān)鍵是如何將信息融合模塊、數(shù)據(jù)同化模塊、源項(xiàng)反演模塊、智能尋優(yōu)算法設(shè)計(jì)并編譯至求解器中。因此，通過關(guān)鍵算法設(shè)計(jì)、基礎(chǔ)參數(shù)和數(shù)據(jù)線索表的制備與讀取、算法庫和求解器的編譯3個(gè)步驟，實(shí)現(xiàn)ChdpFOAM求解器的開發(fā)。

1.2 關(guān)鍵算法設(shè)計(jì)

1.2.1 信息融合算法設(shè)計(jì)

信息融合算法用于對外部監(jiān)測數(shù)據(jù)進(jìn)行處理和融合分析。針對外場環(huán)境面臨化學(xué)危害濃度的差異性大、溫濕度條件不可控等問題，采用聯(lián)合分布函數(shù)對化學(xué)危害濃度和濕度進(jìn)行漏警率影響分析，構(gòu)建傳感器漏警率與當(dāng)前位置濃度、濕度的關(guān)系模型；根據(jù)環(huán)境中汽油、柴油、煤油等對報(bào)警裝備的干擾影響，構(gòu)建單個(gè)傳感器虛警率模型；根據(jù)傳感器分布規(guī)律，設(shè)計(jì)多個(gè)傳感器分布式信息融合規(guī)則，構(gòu)建基于貝葉斯的信息融合判決規(guī)則。信息融合模塊中的模型算法已在文獻(xiàn)[13-14]中詳細(xì)介紹，最終的信息融合判決規(guī)則為

(1)

式中：k為地理網(wǎng)格坐標(biāo)點(diǎn)的個(gè)數(shù)或監(jiān)測節(jié)點(diǎn)個(gè)數(shù)；ui為單個(gè)監(jiān)測節(jié)點(diǎn)的局部判決；PMi、PFi分別表示監(jiān)測節(jié)點(diǎn)i的漏警和虛警概率，PMi=Pr(ui=0|H1)，PFi=Pr(ui=1|H0)；u0為融合后的判決；η為門限。當(dāng)(1)式的左邊大于等于右邊時(shí)，融合輸出u0=1，即認(rèn)為H1成立，即危害事件成立，否則H0成立，即不存在危害事件。

1.2.2 數(shù)據(jù)同化算法設(shè)計(jì)

數(shù)據(jù)同化算法用于對產(chǎn)生化學(xué)危害濃度分布進(jìn)行修正，其數(shù)據(jù)來源為信息融合模塊(觀測數(shù)據(jù))和PISOFOAM計(jì)算模塊(模擬數(shù)據(jù))。相對于其他同化方法，4DVAR實(shí)現(xiàn)了整個(gè)同化時(shí)間窗口觀測數(shù)據(jù)的同化，預(yù)測精度更高，適用范圍更廣，在解決化學(xué)危害擴(kuò)散預(yù)測中具有明顯的優(yōu)勢。本文采用4DVAR作為數(shù)據(jù)同化算法，構(gòu)建適用于化學(xué)危害擴(kuò)散的目標(biāo)泛函梯度模型、數(shù)值預(yù)報(bào)模式及切線性算子?；瘜W(xué)危害擴(kuò)散目標(biāo)泛函模型為

(2)

式中：C為擴(kuò)散空間危害濃度分布；n為危害源數(shù)量；cj,0為第j個(gè)危害源初始時(shí)刻危害物濃度狀態(tài)向量；cb為背景場向量；B為背景誤差協(xié)方差矩陣；T為同化時(shí)間窗口；Ht為t時(shí)刻觀測算子；cj,t為第j個(gè)危害源t時(shí)刻擴(kuò)散濃度狀態(tài)向量；yt為t時(shí)刻危害物濃度觀測向量；R為觀測誤差協(xié)方差矩陣。

(3)

考慮化學(xué)危害擴(kuò)散主要影響因素為對流項(xiàng)，為降低擴(kuò)散預(yù)測誤差，采用泰勒級數(shù)迎風(fēng)差分構(gòu)建各監(jiān)測節(jié)點(diǎn)處的數(shù)值預(yù)報(bào)算子。t時(shí)刻監(jiān)測節(jié)點(diǎn)(x,y,z)處的數(shù)值預(yù)報(bào)算子表達(dá)式為

(4)

根據(jù)t時(shí)刻各個(gè)監(jiān)測節(jié)點(diǎn)的數(shù)值預(yù)報(bào)算子，即可得到t時(shí)刻整個(gè)狀態(tài)向量的數(shù)值預(yù)報(bào)模式Mt=(M1,t,M2,t,…,Mk,t)。同理可得不同時(shí)刻的數(shù)值預(yù)報(bào)模式。假設(shè)t-1時(shí)刻各個(gè)監(jiān)測節(jié)點(diǎn)的危害物預(yù)測濃度依次為ct-1=(c1,t-1,c2,t-1,…,ck,t-1)，根據(jù)切線性算子定義，得到t時(shí)刻數(shù)值預(yù)報(bào)切線性算子Dt：

(5)

1.2.3 源項(xiàng)反演算法設(shè)計(jì)

源項(xiàng)反演算法用于對真實(shí)源強(qiáng)與爆炸位置進(jìn)行反向推演。通過以下3個(gè)步驟構(gòu)建源項(xiàng)反演算法：1)根據(jù)等值線(面)梯度快速估算生成危害源位置的初始猜測，如圖2所示；2)將初始猜測作為初始源項(xiàng)信息進(jìn)行擴(kuò)散模擬仿真，綜合場景中不同位置布設(shè)的觀測節(jié)點(diǎn)陣列的多組觀測值進(jìn)行數(shù)據(jù)同化；3)將同化后的預(yù)測數(shù)據(jù)與觀測數(shù)據(jù)進(jìn)行相對均方根誤差(rRMSE)[15-16]計(jì)算分析，通過智能尋優(yōu)算法實(shí)現(xiàn)循環(huán)迭代，最終獲得rRMSE極小下的危害源位置和強(qiáng)度。圖2中，A為梯度下降基準(zhǔn)點(diǎn)，C、D、E為梯度下降中間節(jié)點(diǎn)，F(xiàn)為危害源位置初始猜測，B為A點(diǎn)等值線的極點(diǎn)。rRMSE計(jì)算模型為

(6)

圖2 追蹤格網(wǎng)數(shù)據(jù)估算等值線的梯度方向Fig.2 Gradient direction of contour line estimatedfrom tracking grid data

1.2.4 智能尋優(yōu)算法設(shè)計(jì)

智能尋優(yōu)算法用于對數(shù)據(jù)同化算法和源項(xiàng)反演算法進(jìn)行快速求解最優(yōu)值，以求得極小目標(biāo)泛函梯度和最優(yōu)源項(xiàng)信息。改進(jìn)花朵授粉算法(IFPA)是Zhu等[17]提出的一種高效智能尋優(yōu)算法，該算法將佳點(diǎn)集種群初始化、Deb可行性比較法、ε約束法融入花朵授粉算法(FPA)，全面提高了FPA的全局尋優(yōu)能力和處理不同約束問題的靈活性。與此同時(shí)，相比于粒子群優(yōu)化(PSO)算法、人工蜂群(ABC)算法、協(xié)同進(jìn)化差分進(jìn)化算法(CDEA)、遺傳算法(GA)等方法，IFPA不受可行域不連通的限制，在求解過程中不需要目標(biāo)函數(shù)的梯度或?qū)?shù)信息，為解決高維、多極值、不連續(xù)、不可微等復(fù)雜系統(tǒng)約束優(yōu)化問題提供了新的思路和手段。因此，本文選取IFPA作為智能尋優(yōu)算法。

1.3 求解器的編譯

ChdpFOAM求解器的編譯分為以下8個(gè)步驟：1)復(fù)制目錄PISOFOAM到新的位置，新目錄名為ChdpFOAM；2)將PISOFOAM.C改名為ChdpFOAM.C；3)刪除依賴文件PISOFOAM.dep；4)修改編譯文件files，進(jìn)入Make目錄，打開files文件，將EXE=$(FOAM_APPBIN)/ PISOFOAM改為EXE=$(FOAM_USER_APPBIN)/ChdpFOAM；5)刪除原來的obj文件；6)在creatFiles.H文件中建立濃度場預(yù)測模型、信息融合算法、數(shù)據(jù)同化算法、源項(xiàng)反演算法、智能尋優(yōu)算法；7)在ChdpFOAM.C中，通過fvScalarMatrix、fvVectorMatrix和.solve()工具加入濃度微分方程(CHEqn)、信息融合(u0Eqn)算法、數(shù)據(jù)同化(J4DEqn)算法、源項(xiàng)反演(rRMSEEqn)算法、智能尋優(yōu)(CSEqn、QEqn)算法；8)編譯wmake，實(shí)現(xiàn)算法模型間數(shù)據(jù)格式和接口設(shè)置。編譯后的ChdpFOAM求解器計(jì)算過程如圖3所示。

圖3 ChdpFOAM求解器計(jì)算過程Fig.3 Calculation process of ChdpFOAM solver

2 ChdpFOAM求解器的適用性驗(yàn)證

以2018年福州大學(xué)某次六氟化硫(SF6)示蹤實(shí)驗(yàn)為例，通過構(gòu)建虛擬地理三維場景，對釋放源信息、氣象信息、擴(kuò)散空間等進(jìn)行設(shè)置或剖分，通過ChdpFOAM求解器對擴(kuò)散過程進(jìn)行模擬和求解。

2.1 SF6示蹤實(shí)驗(yàn)基本情況介紹

2018年11月8日上午6時(shí)，在福州大學(xué)進(jìn)行了大氣擴(kuò)散示蹤實(shí)驗(yàn)，施放示蹤物為SF6，風(fēng)向?yàn)槲鞅憋L(fēng)，風(fēng)速達(dá)中性穩(wěn)定度，釋放量為4.65 kg，釋放時(shí)間為80 s。本次實(shí)驗(yàn)場地空間約為1.0 km×1.5 km，共設(shè)立2個(gè)固定氣象站、7個(gè)便攜式氣象觀測點(diǎn)和120個(gè)不規(guī)則取樣點(diǎn)，以實(shí)時(shí)觀測風(fēng)速、溫度、濕度、氣壓和示蹤劑濃度等數(shù)據(jù)，釋放點(diǎn)和監(jiān)測節(jié)點(diǎn)布置如圖4所示。

圖4 實(shí)驗(yàn)點(diǎn)位布置Fig.4 Layout of experimental points

2.2 實(shí)驗(yàn)?zāi)M與分析

2.2.1 實(shí)驗(yàn)參數(shù)設(shè)置

采用ChdpFOAM求解器對整個(gè)擴(kuò)散過程進(jìn)行模擬計(jì)算，空氣設(shè)為主相，SF6作為第2相，各參數(shù)設(shè)置為：SF6密度6.52 kg/m3，動(dòng)力黏度0.000 014 2 Pa·s。智能尋優(yōu)算法中參數(shù)設(shè)置如下：種群規(guī)模N=200，授粉方式轉(zhuǎn)換概率P=0.8，最大迭代次數(shù)T=500。

2.2.2 網(wǎng)格剖分及邊界條件設(shè)定

通過blockMesh剖分、snappyHexMesh剖分和正交式非均勻交錯(cuò)網(wǎng)格剖分，對擴(kuò)散空間中的釋放源、下墊面和障礙物附近進(jìn)行加密網(wǎng)格劃分，如圖5所示。圖5中，擴(kuò)散空間由x軸正向1 500 m、y軸正向1 000 m、z軸正向100 m的截取方式創(chuàng)建，網(wǎng)格分辨率設(shè)為10 m×10 m。

圖5 擴(kuò)散空間網(wǎng)格剖分Fig.5 Mesh generation of diffusion space

在網(wǎng)格剖分的基礎(chǔ)上，對危害物的擴(kuò)散初始邊界條件進(jìn)行設(shè)定：根據(jù)各氣象觀測點(diǎn)數(shù)據(jù)，構(gòu)建擴(kuò)散空間內(nèi)初始風(fēng)場信息，如圖6所示。危害物釋放點(diǎn)坐標(biāo)為(-527 m，597 m)，釋放高度為10 m，釋放強(qiáng)度為58 g/s，釋放持續(xù)時(shí)間為80 s；入口湍流強(qiáng)度為4.02%，入口面積為10 m×10 m，出口湍流強(qiáng)度為4.15%，外部氣壓為101 920 Pa，濕度24%。初始邊界條件設(shè)定如表1所示。

表1 初始邊界條件設(shè)定Tab.1 Initial boundary condition setting

2.2.3 實(shí)驗(yàn)過程模擬

在實(shí)驗(yàn)參數(shù)設(shè)置和計(jì)算仿真的基礎(chǔ)上，對示蹤劑在不同時(shí)刻、不同高度、不同剖面的濃度分布情況進(jìn)行模擬和可視化展示。

2.2.3.1 不同時(shí)刻擴(kuò)散過程模擬

圖7所示為沿風(fēng)向切面10 s、30 s、1 min、2 min、3 min、4 min時(shí)的示蹤劑擴(kuò)散濃度分布模擬狀況。

圖7 沿風(fēng)向方向各時(shí)刻示蹤劑擴(kuò)散濃度分布效果Fig.7 Tracer diffusion concentration at each time along the wind direction

由圖7可知，隨著時(shí)間的增加，示蹤劑沿風(fēng)向方向逐漸蔓延至整個(gè)擴(kuò)散空間；在遇到建筑物等干擾后，示蹤劑會(huì)從兩側(cè)或上方隨氣流輸送；在建筑物密集區(qū)域，示蹤劑的擴(kuò)散速度會(huì)減緩，部分地帶出現(xiàn)滯留現(xiàn)象。

2.2.3.2 不同高度擴(kuò)散模擬分析

以t=4 min為例，對示蹤劑在不同高度h的空間分布狀況進(jìn)行模擬和分析，如圖8所示。

圖8 不同高度下示蹤劑擴(kuò)散濃度分布效果Fig.8 Tracer diffusion concentration distribution at different heights

由圖8可知，2 m高度時(shí)，示蹤劑擴(kuò)散受地形、建筑物和樹木的影響非常明顯，部分區(qū)域示蹤劑存在明顯滯留現(xiàn)象，危害區(qū)域標(biāo)識呈雜亂無序特征；隨著高度的增加，地形和建筑物對示蹤劑擴(kuò)散的影響逐漸減小，危害區(qū)域也逐漸呈規(guī)則化特征。表明小尺度空間內(nèi)危害物的擴(kuò)散受下墊面影響明顯，采用計(jì)算流體力學(xué)開展化學(xué)危害擴(kuò)散預(yù)測優(yōu)勢十分突出。

2.2.3.3 不同剖面擴(kuò)散模擬分析

以t=4 min為例，對示蹤劑在平行于風(fēng)向、平行于地面、垂直于風(fēng)向等3種剖面下的空間分布狀況進(jìn)行模擬和分析，如圖9所示。

圖9 不同剖面下示蹤劑擴(kuò)散濃度分布效果Fig.9 Tracer diffusion concentration distribution in different sections

圖10 4DVAR同化空間軌跡圖Fig.10 4DVAR assimilation space trajectory

由圖9(a)可知，沿著風(fēng)向方向，隨著擴(kuò)散區(qū)域的增大，除部分建筑物密集區(qū)域外，其余區(qū)域示蹤劑濃度呈逐漸降低態(tài)勢。由圖9(b)可知，同一高度下，示蹤劑濃度分布主要受風(fēng)向和下墊面作用。由圖9(c)可知，同一區(qū)域內(nèi)，不同高度的示蹤劑濃度分布也相差較大，示蹤劑的擴(kuò)散不僅受下墊面影響，受風(fēng)向風(fēng)速和湍流的作用也十分明顯。綜合3個(gè)剖面可知，化學(xué)危害擴(kuò)散預(yù)測不僅要按地域劃分，更要按空間分層，對不同維度危害態(tài)勢進(jìn)行標(biāo)注和顯示。

2.2.4 ChdpFOAM求解器計(jì)算仿真

在實(shí)驗(yàn)參數(shù)設(shè)置基礎(chǔ)上，按照ChdpFOAM求解器計(jì)算流程，對示蹤劑擴(kuò)散實(shí)驗(yàn)進(jìn)行模擬仿真，求解經(jīng)過數(shù)據(jù)同化和智能尋優(yōu)后的流場和濃度場(源項(xiàng)已知)。由于釋放量較小，持續(xù)時(shí)間較短，僅將4 min以內(nèi)觀測數(shù)據(jù)作為校驗(yàn)對象。圖10所示為4DVAR同化空間軌跡圖，圖11所示為IFPA尋優(yōu)變化曲線。

圖11 IFPA尋優(yōu)變化曲線Fig.11 IFPA optimization curve

由圖10可知，4DVAR同化窗口基本覆蓋整個(gè)搜索空間，可通過多組觀測濃度場數(shù)據(jù)對預(yù)測濃度場數(shù)據(jù)進(jìn)行同化處理，達(dá)到求解器同化目的；但是，受實(shí)驗(yàn)限制，觀測濃度場數(shù)據(jù)較少，同化計(jì)算與真實(shí)環(huán)境下數(shù)據(jù)量和計(jì)算量均相差較大，算法的穩(wěn)定性還需進(jìn)一步進(jìn)行驗(yàn)證。

由圖11可知，5次尋優(yōu)計(jì)算在迭代200～300次時(shí)均基本達(dá)到穩(wěn)定收斂狀態(tài)，J4DVAR(C)可逐漸降低至10%左右，說明IFPA的整體收斂性較強(qiáng)，預(yù)測結(jié)果基本滿足現(xiàn)實(shí)需求，可通過減少迭代次數(shù)提高算法的計(jì)算速率。

2.3 實(shí)驗(yàn)結(jié)果分析

2.3.1 誤差原因分析

由模擬實(shí)驗(yàn)可知，通過ChdpFOAM求解器可使J4DVAR(C)值降至10%左右，誤差原因主要以下3個(gè)方面：

1)尋優(yōu)算法計(jì)算誤差。IFPA只能通過循環(huán)迭代求得相對最優(yōu)值，很難類似窮舉法得到最優(yōu)值。此類誤差可通過多次獨(dú)立仿真進(jìn)行降低，但是不可能避免的。通過對上述實(shí)驗(yàn)進(jìn)行20次獨(dú)立實(shí)驗(yàn)，最終J4DVAR(C)最小值為9.54%。

2)觀測儀器測量誤差。示蹤劑取樣設(shè)備、氣象測量設(shè)備、示蹤劑分析設(shè)備等測量偏差也是導(dǎo)致計(jì)算誤差產(chǎn)生的原因。為此，對0.522 μg/m3濃度的SF6進(jìn)行了7次分析，相對標(biāo)準(zhǔn)偏差為4.1%，如表2所示。

表2 SF6分析相對標(biāo)準(zhǔn)偏差Tab.2 Relative standard deviation of SF6 analysis

3)背景場誤差。造成實(shí)驗(yàn)預(yù)測模擬誤差的另一個(gè)重要原因是背景場誤差，主要由下墊面干擾誤差和風(fēng)場誤差構(gòu)成。下墊面干擾誤差與地理模型精細(xì)程度、氣流與建筑物的干涉等息息相關(guān)。風(fēng)場誤差主要由監(jiān)測節(jié)點(diǎn)的數(shù)量和布局不夠系統(tǒng)均勻造成。

2.3.2 ChdpFOAM求解器優(yōu)劣性分析

為驗(yàn)證ChdpFOAM求解器優(yōu)劣性，按照均勻分布、覆蓋面全的規(guī)則選取10個(gè)監(jiān)測節(jié)點(diǎn)，以90 s、120 s、180 s、240 s 4個(gè)時(shí)刻數(shù)據(jù)作為分析數(shù)據(jù)，分別通過ChdpFOAM求解器和PISOFOAM求解器計(jì)算求解，分析二者與真實(shí)觀測值之間的差異性，結(jié)果如表3所示。各監(jiān)測節(jié)點(diǎn)坐標(biāo)依次為：(-420 m，570 m)、(-170 m，600 m)、(-190 m，820 m)、(-340 m，410 m)、(-320 m，100 m)、(-400 m，450 m)、(-330 m，260 m)、(-360 m，-60 m)、(-330 m，130 m)、(-350 m，-210 m)。濃度差異計(jì)算表達(dá)式為

(7)

式中：d1為ChdpFOAM求解器模擬預(yù)測濃度值與真實(shí)觀測濃度值的差異；cC為ChdpFOAM求解器模擬預(yù)測濃度值；cT為真實(shí)觀測濃度值；d2為PISOFOAM求解器模擬預(yù)測濃度值與真實(shí)觀測濃度值的差異；cp為PISOFOAM求解器模擬預(yù)測濃度值；d3為ChdpFOAM求解器相對PISOFOAM求解器模擬精度提高程度。通過數(shù)據(jù)分析，將ChdpFOAM、PISOFOAM求解器擴(kuò)散模擬與真實(shí)觀測結(jié)果進(jìn)行對比分析，如圖12和圖13所示。

由表3和圖12、圖13可知：d1最小值(4.71%)出現(xiàn)在6號監(jiān)測節(jié)點(diǎn)(120 s)，最大值(34.73%)出現(xiàn)在8號監(jiān)測節(jié)點(diǎn)(240 s)，平均誤差百分比為12.98%；d2最小值(7.76%)出現(xiàn)在10號監(jiān)測節(jié)點(diǎn)(180 s)，最大值(33.98%)出現(xiàn)在5號監(jiān)測節(jié)點(diǎn)(180 s)，平均誤差百分比為18.95%；d3最小值(-0.20%)和最大值(16.43%)均出現(xiàn)在240 s時(shí)刻(1號、3號監(jiān)測節(jié)點(diǎn))，平均誤差差異為5.97%。從模擬結(jié)果角度分析，在下墊面相對平坦、建筑物少的區(qū)域(6號監(jiān)測節(jié)點(diǎn)、10號監(jiān)測節(jié)點(diǎn))，ChdpFOAM和PISOFOAM求解器計(jì)算誤差均較低；在距建筑物相對密集的區(qū)域(8號監(jiān)測節(jié)點(diǎn)、5號監(jiān)測節(jié)點(diǎn))，ChdpFOAM和PISOFOAM求解器計(jì)算誤差均較高；相比于PISOFOAM求解器，ChdpFOAM求解器的計(jì)算精度總體提升6%左右；180 s內(nèi)ChdpFOAM求解器計(jì)算精度明顯高于PISOFOAM求解器，但是240 s時(shí)刻同化效果比較微弱。

表3 ChdpFOAM、PISOFOAM求解器擴(kuò)散模擬與真實(shí)觀測結(jié)果對比Tab.3 Comparison of the simulated diffusion results ChdpFOAM and PISOFOAM solvers and the real observed results

圖12 ChdpFOAM、PISOFOAM求解器擴(kuò)散模擬與真實(shí)觀測數(shù)據(jù)分析Fig.12 Analysis of the simulated diffusion data of ChdpFOAM and PISOFOAM solvers and the true observed data

圖13 ChdpFOAM和PISOFOAM求解器預(yù)測誤差分析Fig.13 Prediction error analysis between ChdpFOAM and PISOFOAM

為進(jìn)一步分析植入模塊后對求解器計(jì)算效率的影響，分別對不同迭代次數(shù)、不同時(shí)間步長下ChdpFOAM和PISOFOAM求解器需要的計(jì)算時(shí)間進(jìn)行統(tǒng)計(jì)分析，如表4所示，其中d4表示ChdpFOAM和PISOFOAM求解器的耗時(shí)比。

表4 ChdpFOAM和PISOFOAM求解器不同步長模擬時(shí)間消耗統(tǒng)計(jì)Tab.4 Time consumption statistics of ChdpFOAM andPISOFOAM solvers at different step sizes

由表4可見：以模擬240 s擴(kuò)散為前提，使用ChdpFOAM和PISOFOAM求解器對不同時(shí)間步長時(shí)的擴(kuò)散進(jìn)行模擬，時(shí)間步長越短、迭代次數(shù)越多，二者的耗時(shí)越長；以30 s步長為例，采用PISOFOAM求解器耗時(shí)729 s，采用ChdpFOAM求解器耗時(shí)967 s， PISOFOAM求解器的計(jì)算效率相對更高；從耗時(shí)比角度分析，ChdpFOAM求解器與PISOFOAM求解器的平均耗時(shí)比為1.33，計(jì)算效率下降了約32.57%。說明信息融合、數(shù)據(jù)同化、智能優(yōu)化等模塊的運(yùn)行，未對求解器計(jì)算效率造成明顯下降，基本滿足實(shí)際需求。

3 結(jié)論

本文針對化學(xué)危害擴(kuò)散受地形、障礙物、湍流等影響顯著的特點(diǎn)，基于OpenFOAM軟件開發(fā)平臺(tái)設(shè)計(jì)適用于化學(xué)危害擴(kuò)散預(yù)測的ChdpFOAM求解器，并對其適用性進(jìn)行實(shí)驗(yàn)驗(yàn)證。通過實(shí)驗(yàn)驗(yàn)證分析，得到如下主要結(jié)論：

1)小尺度空間內(nèi)，危害物的擴(kuò)散受下墊面作用明顯，ChdpFOAM求解器具有較強(qiáng)的實(shí)用性。在下墊面相對平坦、建筑物少的區(qū)域，干擾因素較少，ChdpFOAM求解器的預(yù)測誤差相對較低；在建筑物密集區(qū)域，SF6的擴(kuò)散速度會(huì)減緩，部分地帶出現(xiàn)滯留現(xiàn)象，ChdpFOAM求解器的求解難度大大增加，可能導(dǎo)致預(yù)測精度的降低。

2)沿風(fēng)向方向，由于SF6的密度約為空氣的5倍，易受沉降效應(yīng)聚集形成近地面毒劑云團(tuán)；遠(yuǎn)離釋放點(diǎn)位置，經(jīng)過建筑物和下墊面的湍流擾動(dòng)，部分SF6會(huì)發(fā)生抬升導(dǎo)致濃度下降，經(jīng)過建筑物后又逐漸向近地面沉降。

3)相對于PISOFOAM求解器，ChdpFOAM求解器在計(jì)算精度約提升6%，求解精度改進(jìn)效果明顯。

4)相對于PISOFOAM求解器，ChdpFOAM求解器計(jì)算效率下降約32.57%，下一步可通過模塊優(yōu)化和內(nèi)部接口設(shè)計(jì)提升計(jì)算效率。