基于正交試驗(yàn)法的高斯擴(kuò)散模型預(yù)測精度研究*

馬寒簫,邢志祥,吳 凡,吳 潔

(常州大學(xué) 安全科學(xué)與工程學(xué)院,江蘇 常州 213164)

0 引言

近年來,我國化工行業(yè)迅速崛起,化工園區(qū)大量涌現(xiàn),為地方經(jīng)濟(jì)帶來貢獻(xiàn)的同時(shí)也伴隨一定安全隱患。大型化工廠氣體泄漏引發(fā)的事故頻發(fā)[1-2],尤其是化工園區(qū)內(nèi)有毒有害氣體泄漏導(dǎo)致的中毒、爆炸事故造成嚴(yán)重的人員傷亡和財(cái)產(chǎn)損失。為保障園區(qū)內(nèi)人員安全,研究開放空間下有毒有害氣體在大氣中的濃度分布,有效地預(yù)測有毒有害物質(zhì)的擴(kuò)散規(guī)律,可為后續(xù)應(yīng)急救援工作的順利開展奠定基礎(chǔ)。

高斯擴(kuò)散模型作為研究輕氣擴(kuò)散規(guī)律的經(jīng)驗(yàn)?zāi)Ｐ?具有處理模糊性和隨機(jī)性的能力,是研究定性概念和定量數(shù)值轉(zhuǎn)換的數(shù)學(xué)模型[3]。目前,國內(nèi)外學(xué)者對于高斯模型的預(yù)測精度研究較少,大部分研究集中于應(yīng)用模型進(jìn)行事故場景分析。例如,倪健等[4]考慮溫度、風(fēng)速、季風(fēng)風(fēng)向變化的影響,采用高斯擴(kuò)散模型分析污染源的數(shù)量、濃度大小以及位置分布關(guān)系;熊立春等[5]針對現(xiàn)有高斯模型進(jìn)行液氨泄漏擴(kuò)散模擬時(shí)未考慮氣體空間疊加效應(yīng)的問題,提出引入時(shí)間疊加因素的高斯擴(kuò)散模型;Shi等[6]運(yùn)用高斯擴(kuò)散模型對液氫泄漏擴(kuò)散特性進(jìn)行研究;Marro等[7]采用積分煙羽上升模型預(yù)測煙羽的中心軌跡;Pournazeri等[8]修改守恒方程,以考慮整體羽流上升模型中的上升氣流和下降氣流;Tao等[9]研究兩相羽流動量主導(dǎo)型和浮力主導(dǎo)型區(qū)域特征,利用數(shù)值模擬結(jié)果修改積分煙羽上升模型的系數(shù);Liu等[10]認(rèn)為有必要對簡單模型進(jìn)行修改,以適應(yīng)更復(fù)雜的環(huán)境。由于大氣穩(wěn)定度、氣體種類、氣體初始密度這3個(gè)因素對氣體擴(kuò)散濃度分布影響較大,因此,研究這些參數(shù)對高斯模型預(yù)測精度的影響具有實(shí)際意義。

因此,本文擬通過正交試驗(yàn)法研究大氣穩(wěn)定度、氣體種類、氣體初始密度這3個(gè)因素對高斯模型預(yù)測精度的影響,通過將FDS軟件獲得的模擬值與MATLAB軟件計(jì)算得到的模型預(yù)測值進(jìn)行對比,引入歸一化均方差NMSE,對比模擬值與預(yù)測值的誤差,以更好地展現(xiàn)高斯模型預(yù)測精度的影響因素和不同影響因素的顯著性,并結(jié)合模擬結(jié)果分析其原因。

1 試驗(yàn)部分

1.1 三因素五水平正交表

正交試驗(yàn)法作為研究多因素多水平的試驗(yàn)設(shè)計(jì)方法,具有均衡分散、數(shù)據(jù)計(jì)算簡單、水平整齊可比的優(yōu)點(diǎn)。本文采用三因素五水平正交表L25(53)進(jìn)行模擬試驗(yàn),如表1所示。

表1 因素水平表Table 1 Levels of factors

1.2 模型建立及邊界條件設(shè)置

高斯氣體擴(kuò)散模型作為較早的經(jīng)驗(yàn)?zāi)Ｐ痛嬖谄湟欢ǖ膬?yōu)勢以及局限性。高斯模型分為高斯煙羽模型和高斯煙團(tuán)模型,當(dāng)輕氣連續(xù)泄漏時(shí),應(yīng)采用高斯煙羽模型,其氣體濃度分布規(guī)律如式(1)所示:

(1)

式中:h表示羽流中心線距離地面的有效高度,m;c(x,y,z)表示質(zhì)量濃度(以下簡稱濃度),kg/m3;σy,σz表示擴(kuò)散系數(shù),m;z表示垂直高度,m;q表示質(zhì)量流量,kg/s;Ua表示風(fēng)速,m/s。本文所述濃度均為質(zhì)量濃度。

國內(nèi)外對于高斯模型羽流中心線有效高度的計(jì)算公式差異較大,因此為確定其最優(yōu)高度,將模擬值與預(yù)測值濃度變化曲線進(jìn)行對比,選取二者濃度變化曲線峰值相等時(shí)的羽流中心線高度為最優(yōu)高度,如圖1所示。采用FDS軟件建立的數(shù)值模型長500 m,寬300 m,高50 m,泄漏源高度設(shè)置為4 m,泄漏面積2 m×2 m,泄漏量20 kg/s,位置距左側(cè)邊界120 m,距右側(cè)邊界380 m。在泄漏源的中心線距離地面10 m高度處設(shè)置濃度測點(diǎn)用于記錄泄漏氣體10 min內(nèi)的平均濃度。

圖1 預(yù)測值與模擬值濃度對比Fig.1 Comparison of theoretical and simulated concentrations

圖2 FDS模擬模型Fig.2 FDS simulation model

1.3 模型驗(yàn)證

采用試驗(yàn)方法研究有害氣體和煙氣在大空間的擴(kuò)散具有成本高、誤差大等弊端,因此,本文采用數(shù)值模擬對各類氣體在大空間的擴(kuò)散開展研究。為確保數(shù)值模擬的準(zhǔn)確性,模擬前將各大氣穩(wěn)定度下風(fēng)速與溫度分布模擬結(jié)果與相關(guān)試驗(yàn)數(shù)據(jù)進(jìn)行對比,大氣穩(wěn)定度C、F等級下風(fēng)速、溫度隨高度變化的對比結(jié)果如圖3所示。試驗(yàn)數(shù)據(jù)來源于2006年美國運(yùn)輸部管道和危險(xiǎn)材料安全管理局下屬的桑迪亞國家試驗(yàn)室開展的大氣擴(kuò)散試驗(yàn)[11-12]。

圖3 C、F等級大氣穩(wěn)定度下風(fēng)速及溫度隨高度的變化Fig.3 Comparison of temperature and wind speed under C and F levels

由圖3可知,采用FDS軟件模擬得到的溫度、風(fēng)速變化與試驗(yàn)值吻合較好,說明FDS軟件對于不同大氣穩(wěn)定度下的風(fēng)速分布以及溫度分布具有較高的預(yù)測精度,可用于預(yù)測泄漏氣體在不同大氣穩(wěn)定度下的擴(kuò)散分布。

2 模擬試驗(yàn)結(jié)果及分析

2.1 試驗(yàn)結(jié)果

為研究高斯擴(kuò)散模型預(yù)測精度的主要影響因素,本文采用正交試驗(yàn)的方法選取B、C、D、E、F 5種大氣穩(wěn)定度以及H2、CH4、CO、SO2、C7H85種化工園區(qū)中常見氣體,并分別設(shè)置0.6 kg/m3、0.7 kg/m3、0.8 kg/m3、0.9 kg/m3、1.0 kg/m35種常見氣體初始密度。為了更好的對不同工況中模擬值與理論預(yù)測值進(jìn)行對比,引入歸一化均方誤差NMSE來量化模擬值與預(yù)測值的吻合程度,如式(2)所示:

(2)

式中:C0表示濃度預(yù)測值,g/m3;CP表示濃度模擬值,g/m3。

NMSE值越小,表示模擬值與預(yù)測值越接近,NMSE值越大,表示模擬值與預(yù)測值相差越大。正交試驗(yàn)結(jié)果如表2所示。

2.2 極差分析及方差分析

對模擬后的25組工況進(jìn)行極差分析,分別計(jì)算每個(gè)因素的極差,結(jié)果見表3。

其中,ki表示NMSE值的求和平均值;R為極差,其值越大,表明該因素對指標(biāo)的影響較大,一般為主要因素;反之,表明該因素對指標(biāo)的影響較小,一般為次要因素[13]。極差計(jì)算如式(3)所示:

圖4為硅鋁比為3.0、pH值為9.73、導(dǎo)向劑用量為3%、老化溫度為60 ℃、晶化溫度100 ℃、晶化時(shí)間為24 h時(shí)，合成的13X分子篩的粒度分布圖。13X分子篩的粒徑主要由兩部分組成：<1 μm和1～10 μm；顆粒體積占顆粒總體積的50%時(shí)，分子篩的粒徑為2.8 μm。

R=kmax-kmin

(3)

極差分析相較于方差分析較簡單且計(jì)算量較小,但極差分析存在一定局限性。極差法沒有將各因素的水平變化與試驗(yàn)誤差所引起的數(shù)據(jù)波動嚴(yán)格區(qū)分開,方差分析則可以通過計(jì)算用數(shù)字量化各因素和誤差的影響[14-15]。因此,本文將進(jìn)一步采用方差分析進(jìn)行計(jì)算,以確保結(jié)論的準(zhǔn)確性,如式(4)～(6)所示:

f=m-1

(4)

(5)

(6)

式中:m為水平數(shù),取值5;n為試驗(yàn)次數(shù),取值25;f為自由度,取值4;sj為離差平方和;E為均方。

方差分析結(jié)果如表4所示。

表4 方差分析Table 4 Analysis of variance

由表3和表4可知,氣體泄漏擴(kuò)散濃度的影響因素中對結(jié)果影響顯著性從大到小依次為大氣穩(wěn)定度>氣體種類>氣體初始密度。

2.3 氣體擴(kuò)散濃度變化分析

為進(jìn)一步研究大氣穩(wěn)定度、氣體種類、氣體初始密度對氣體濃度擴(kuò)散的影響,從25組工況中篩選典型的4組工況進(jìn)行分析?；谇拔腘MSE值,本文選取大氣穩(wěn)定度、氣體種類、氣體初始密度中吻合最好、吻合較好、吻合最差3種工況進(jìn)行二維云圖濃度分析,并采用Tecplot軟件進(jìn)行對比分析。

2.3.1 大氣穩(wěn)定度

正交試驗(yàn)結(jié)果分析中,只考慮大氣穩(wěn)定度這一因素的影響,氣體擴(kuò)散濃度變化預(yù)測值與模擬值吻合最好的是大氣穩(wěn)定度C等級,吻合最差的是大氣穩(wěn)定度F等級。本文選取第7組工況與第25組工況進(jìn)行分析,如圖4所示。

圖4 第7、25組工況模擬及預(yù)測濃度擴(kuò)散云圖Fig.4 Cloud diagram of simulated and predicted concentration diffusion under conditions of group 7 and group 25

由圖4(a)可知,第7組工況預(yù)測結(jié)果與模擬結(jié)果較為接近,濃度預(yù)測值峰值出現(xiàn)在X軸方向約10 m處,濃度峰值模擬值出現(xiàn)在約15 m處。在高濃度范圍內(nèi),Y軸方向上預(yù)測結(jié)果相較于模擬值濃度偏低,中、低濃度變化趨勢吻合較好。濃度預(yù)測結(jié)果變化梯度10～20 g/m3對應(yīng)X軸40～60 m處,同等濃度變化梯度模擬結(jié)果位于X軸方向40～50 m。采用高斯擴(kuò)散模型預(yù)測的濃度變化梯度下降較快,但其總體預(yù)測范圍與模擬值較為接近,所以在實(shí)際應(yīng)用中仍然能夠發(fā)揮重要作用。

由圖4(b)可知,第25組工況中,在大氣穩(wěn)定度為F等級下,濃度變化趨勢預(yù)測結(jié)果與模擬結(jié)果變化趨勢差別較大。預(yù)測值在X軸方向30 m處擴(kuò)散濃度達(dá)到10 g/m3,70 m處達(dá)到峰值,為70～80 g/m3;模擬值在X軸方向約0 m處氣體擴(kuò)散濃度既達(dá)到10 g/m3,5 m處達(dá)到濃度峰值,為80～90 g/m3。相較于模擬值,預(yù)測值的濃度梯度上升及下降速度較慢。

氣體泄漏過程中處于不同大氣穩(wěn)定度,實(shí)際環(huán)境溫度與風(fēng)速變化分布存在差異,第7組工況是以大氣穩(wěn)定度為C級、初始密度為0.7 kg/m3的CH4氣體為試驗(yàn)對象,在C級大氣穩(wěn)定度下,CH4的初始密度隨環(huán)境溫度與風(fēng)速的變化發(fā)生改變,但該等級下氣體密度的改變與氣體擴(kuò)散的初始密度相差較小。高斯擴(kuò)散模型計(jì)算是取X方向10 m處的風(fēng)速進(jìn)行計(jì)算,在C等級大氣穩(wěn)定度下,10 m及以上的風(fēng)速變化較小,因此該組工況預(yù)測值與模擬值吻合較好。在F等級大氣穩(wěn)定度下,10 m及以上的風(fēng)速呈線性變化,因此,取10 m高度處的風(fēng)速進(jìn)行濃度計(jì)算誤差較大。在F等級大氣穩(wěn)定度下,溫度呈線性變化,這直接影響泄漏氣體的密度變化,導(dǎo)致受環(huán)境影響的氣體密度與氣體初始密度產(chǎn)生較大差異,最終影響模型預(yù)測精度。

此外,即使空氣靜止,大氣中溫差將引起空氣運(yùn)動,穩(wěn)定條件下也存在湍流。在氣體擴(kuò)散過程中,大氣越穩(wěn)定,湍流越容易被抑制,擴(kuò)散能力越弱;反之,大氣越不穩(wěn)定,則熱力湍流發(fā)展旺盛,對流強(qiáng)烈導(dǎo)致擴(kuò)散能力加強(qiáng)。氣體的擴(kuò)散能力直接影響高斯擴(kuò)散模型中的擴(kuò)散系數(shù)選取,不同大氣穩(wěn)定度的擴(kuò)散系數(shù)不同,從而影響最終的濃度變化。

2.3.2 氣體種類及氣體初始密度

本文選取氣體種類及氣體初始密度這2種影響因素下吻合最好以及吻合較好的第12組和第19組工況進(jìn)行分析,如圖5所示。

圖5 第12和第19組工況模擬機(jī)預(yù)測濃度擴(kuò)散云圖Fig.5 Cloud diagram of simulated and predicted concentration diffusion under conditions of group 12 and group 19

由圖5(a)可知,第12組工況濃度變化梯度預(yù)測值相較于模擬值下降慢,且在Y軸方向高濃度區(qū)域預(yù)測值偏低。在D等級大氣穩(wěn)定度下,10 m及以上高度處風(fēng)速變化較小,穩(wěn)定在4.5 m/s～5 m/s,溫度始終保持在24 ℃。模擬設(shè)置CH4氣體初始密度0.6 kg/m3、初始溫度132 ℃,此時(shí)環(huán)境溫度低于初始溫度,擴(kuò)散過程中氣體向環(huán)境放熱,使得密度逐漸增大,更容易受風(fēng)速影響向大氣上層快速擴(kuò)散,導(dǎo)致濃度梯度下降較快。

第19組工況是在E等級大氣穩(wěn)定度下初始密度為0.9 kg/m3的SO2氣體。由圖5(b)可知,濃度預(yù)測值與模擬值差異較大,預(yù)測值初始擴(kuò)散位置在X軸方向明顯滯后預(yù)測值約10 m,預(yù)測值濃度梯度下降慢,高濃度區(qū)域在X軸方向濃度預(yù)測值較低。模擬設(shè)置SO2氣體初始密度為0.9 kg/m3,根據(jù)理想狀態(tài)氣體方程可知此時(shí)溫度為592 ℃。在E等級大氣穩(wěn)定度下,大氣溫度維持在33 ℃左右。當(dāng)氣體泄漏至大氣中時(shí)溫度驟然降低,氣體向周圍環(huán)境放熱,密度逐漸增加至2.5 kg/m3后繼續(xù)以重氣的性質(zhì)進(jìn)行擴(kuò)散,導(dǎo)致氣體無法隨風(fēng)擴(kuò)散更遠(yuǎn)距離,因此,泄漏氣體在X軸方向0 m處迅速擴(kuò)散至最高濃度,濃度梯度受密度變化影響迅速改變。由于氣體擴(kuò)散密度的瞬間改變,導(dǎo)致氣體擴(kuò)散時(shí)產(chǎn)生蒸汽云抬升高度不同,這也直接影響高斯擴(kuò)散模型中的羽流上升高度,從而使高斯擴(kuò)散模型預(yù)測精度降低。由此可以看出,不同氣體種類也會導(dǎo)致氣體重力、密度變化不同,從而影響羽流上升高度并最終影響高斯擴(kuò)散模型的預(yù)測精度。

3 結(jié)論

1)在本文正交試驗(yàn)因素水平范圍內(nèi),高斯擴(kuò)散模型在C等級大氣穩(wěn)定度下氣體擴(kuò)散濃度變化預(yù)測值與模擬值吻合最好,F等級大氣穩(wěn)定度下吻合最差。

2)極差分析和方差分析結(jié)果表明,在采用高斯擴(kuò)散模型預(yù)測氣體濃度變化規(guī)律時(shí),大氣穩(wěn)定度對于氣體濃度擴(kuò)散規(guī)律具有顯著性影響,3種因素對氣體濃度的影響程度為大氣穩(wěn)定度>氣體種類>氣體初始密度。

3)受大氣湍流發(fā)展影響,風(fēng)速及溫度變化曲線差異較大。由于高斯擴(kuò)散模型計(jì)算僅取10 m高度處的風(fēng)速進(jìn)行計(jì)算,因此,風(fēng)速及溫度變化曲線呈線性增加的大氣穩(wěn)定度使得高斯擴(kuò)散模型預(yù)測精度降低。

4)氣體種類不同導(dǎo)致氣體在擴(kuò)散時(shí)的重力與動能有一定差別,從而影響羽流中心線高度的計(jì)算并最終導(dǎo)致模型預(yù)測氣體濃度變化趨勢與模擬結(jié)果差異較大。

5)本文在采用高斯擴(kuò)散模型計(jì)算氣體擴(kuò)散濃度變化規(guī)律時(shí),將模擬值與預(yù)測值峰值相等時(shí)的高度設(shè)為最優(yōu)羽流中心線高度,未來將會在羽流中心線高度方面繼續(xù)深入研究,嘗試找出能夠預(yù)測最優(yōu)羽流中心線高度的預(yù)測模型。