羅振敏 胡騰 李俊 王濤 高有軍 楊勇 王登飛

（1.西安科技大學(xué)安全科學(xué)與工程學(xué)院，西安 710054；2.陜西省工業(yè)過(guò)程安全與應(yīng)急救援工程技術(shù)研究中心，西安 710054；3.陜西省煤火災(zāi)害防治重點(diǎn)實(shí)驗(yàn)室，西安 710054；4.榆林榆川天然氣有限責(zé)任公司，陜西 榆林 719099）

0 引言

近年來(lái)我國(guó)天然氣泄漏著火爆炸事故屢有發(fā)生。2021年10月18日，河北省邯鄲市發(fā)生天然氣泄漏事故，造成3人窒息死亡。2021年10月21日，遼寧省沈陽(yáng)市和平區(qū)一家飯店發(fā)生天然氣泄漏爆炸事故，造成3人死亡、30多人受傷。2021年6月13日，湖北省十堰市發(fā)生重大天然氣爆炸事故，致26人死亡。2021年10月24日，遼寧省瓦房店市一居民樓發(fā)生燃?xì)忾W爆事故，造成2人死亡、7人受傷。因此，研究天然氣在各種工況下的泄漏擴(kuò)散規(guī)律及濃度、速度分布規(guī)律對(duì)人員提前采取防護(hù)措施和制定火災(zāi)爆炸預(yù)防方案具有重要意義。數(shù)值模擬研究具有可操作性強(qiáng)、成本低、應(yīng)用面廣等優(yōu)勢(shì)［1］。

國(guó)內(nèi)外大量學(xué)者對(duì)天然氣泄漏擴(kuò)散進(jìn)行了數(shù)值模擬研究，其中基于CFD模型進(jìn)行研究的有：王春園［2］在FLUENT中對(duì)管道正常情況和泄漏情況分別進(jìn)行了仿真，得到了燃?xì)夤艿佬孤┖蠊艿朗啄﹥啥藲怏w壓力和流量的變化規(guī)律；YUE C J等［3］利用計(jì)算流體動(dòng)力學(xué)軟件FLACS，對(duì)徐州某LNG儲(chǔ)罐泄漏可能造成的多種災(zāi)害進(jìn)行了模擬分析；王文和等［4］以重慶某無(wú)人值守井站為實(shí)際場(chǎng)景，采用ANSYS軟件對(duì)含硫天然氣泄漏擴(kuò)散過(guò)程進(jìn)行模擬?；贏LOHA軟件進(jìn)行研究的有：彭琳等［5］為提高地面天然氣管道泄漏擴(kuò)散范圍預(yù)測(cè)的精度，基于ALOHA事故后果模擬分析和多元回歸預(yù)測(cè)方法建立地面天然氣泄漏擴(kuò)散范圍預(yù)測(cè)模型?；诟咚篃熡鹉Ｐ?、內(nèi)嵌UDM模型的DNV（挪威船級(jí)社）PHAST、SAFETY軟件進(jìn)行研究的有：何寧輝等［6］通過(guò)DNV公司的SAFETI軟件對(duì)于某儲(chǔ)配廠的天然氣罐區(qū)進(jìn)行了泄漏后果的模擬和定量風(fēng)險(xiǎn)分析?；陂_源CFD軟件Open-FOAM進(jìn)行研究的有：周理［7］利用開源CFD計(jì)算程序包OpenFOAM對(duì)激波管內(nèi)的高壓燃?xì)庑孤┳匀棘F(xiàn)象進(jìn)行了數(shù)值模擬。

以上研究都未對(duì)氣云速度矢量在不同方向的分量分布規(guī)律進(jìn)行分析。本文基于FLACS軟件，首先對(duì)文獻(xiàn)［7-8］中的物理實(shí)驗(yàn)得出的天然氣泄漏擴(kuò)散濃度與速度結(jié)果進(jìn)行數(shù)值模擬可靠性驗(yàn)證，然后對(duì)某天然氣儲(chǔ)罐區(qū)進(jìn)行泄漏擴(kuò)散模擬，得出天然氣在一定工況下泄漏擴(kuò)散的濃度、速度矢量分量、等效可燃?xì)庠企w積Q9變化規(guī)律。以期對(duì)處置天然氣儲(chǔ)配氣站場(chǎng)泄漏擴(kuò)散事故具有指導(dǎo)意義。

1 模型建立與場(chǎng)景設(shè)置

本文模擬對(duì)象為某天然氣儲(chǔ)配站，如圖1所示，罐區(qū)共有8個(gè)球形鋼制儲(chǔ)罐（編號(hào)L1—L8）和3間控制室和廠房，每個(gè)球罐體積約為2 000 m3（球罐直徑約16 m，底部距地面5.0 m），球罐與球罐之間布置有扶梯和管道。

圖1 天然氣儲(chǔ)罐區(qū)CFD幾何模型

在FLACS軟件CASD建立相應(yīng)的天然氣儲(chǔ)罐區(qū)CFD幾何模型。假設(shè)正北方向?yàn)?Y，計(jì)算區(qū)域坐標(biāo)為：X（-5～248 m）；Y（-20～115m）；Z（0～35 m）；圖左下角為坐標(biāo)原點(diǎn)。對(duì)CFD模型劃分網(wǎng)格，并對(duì)泄漏源附近網(wǎng)格細(xì)化，在主要研究區(qū)域之外，對(duì)網(wǎng)格拉伸，在細(xì)網(wǎng)格和粗網(wǎng)格之間進(jìn)行平滑操作，最終劃分網(wǎng)格982 688個(gè)。然后對(duì)網(wǎng)格模型進(jìn)行孔隙度計(jì)算。

天然氣在泄漏擴(kuò)散過(guò)程中遵循質(zhì)量守恒方程、能量守恒方程、動(dòng)量守恒方程。

2 濃度與速度結(jié)果可靠性驗(yàn)證

本文分別對(duì)文獻(xiàn)［7-8］中的小尺寸物理實(shí)驗(yàn)進(jìn)行濃度及速度結(jié)果驗(yàn)證。在FLACS中建立與文獻(xiàn)描述完全相同的CFD幾何模型，按文獻(xiàn)所述進(jìn)行泄漏源位置、泄漏方向、泄漏孔徑及其他初始條件的設(shè)置，對(duì)監(jiān)測(cè)點(diǎn)進(jìn)行適當(dāng)刪減，只保留文獻(xiàn)［7］模型中心濃度監(jiān)測(cè)點(diǎn)和文獻(xiàn)［8］4號(hào)速度測(cè)點(diǎn)。濃度結(jié)果及速度結(jié)果的實(shí)驗(yàn)值與模擬值對(duì)比如圖2—圖3所示：

圖2 濃度結(jié)果對(duì)比

圖3 速度結(jié)果對(duì)比

由圖2—圖3可得利用FLACS的模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的變化趨勢(shì)相近，但局部誤差較大，分析誤差原因?yàn)椋涸诂F(xiàn)場(chǎng)實(shí)驗(yàn)過(guò)程中大氣條件發(fā)生一定變化，而數(shù)值模擬過(guò)程中假設(shè)大氣條件不改變。此外，文獻(xiàn)中部分初始條件沒有詳細(xì)列出，可能與模擬中的初始條件有差異，最終導(dǎo)致產(chǎn)生誤差。圖2實(shí)驗(yàn)與模擬數(shù)據(jù)平均相對(duì)誤差為15.36%，圖3實(shí)驗(yàn)與模擬數(shù)據(jù)平均相對(duì)誤差為11.45%，誤差均在20%以內(nèi)，符合工程允許誤差范圍之內(nèi)。綜上所述，F(xiàn)LACS數(shù)值模擬得出的濃度與速度結(jié)果相對(duì)可靠。

3 數(shù)值模擬結(jié)果分析

本文模擬L2儲(chǔ)罐底部發(fā)生泄漏，泄漏源坐標(biāo)為（45 m，50 m，5 m），泄漏孔徑40 mm，泄漏方向?yàn)?Y（正北）。環(huán)境溫度20℃，大氣壓101 325Pa，地面粗糙度0.01，大氣穩(wěn)定度等級(jí)F（由環(huán)境風(fēng)速、云覆蓋面積等天氣情況決定），泄漏時(shí)間為60 s，泄漏停止后持續(xù)擴(kuò)散時(shí)間為30 s，具體工況見表1。

表1 天然氣泄漏擴(kuò)散工況

3.1 風(fēng)速對(duì)氣云擴(kuò)散的影響

圖4—圖7分別為case1、case2、case3、case4條件下泄漏后的天然氣擴(kuò)散達(dá)到穩(wěn)定狀態(tài)時(shí)Z=5 m的XY切面（泄漏源所在平面）天然氣體積分?jǐn)?shù)（1.5%～10%）分布圖、天然氣速度矢量（VVEC）在V方向上的分量（-35～35 m/s）分布圖、天然氣速度矢量（VVEC）在U方向上的分量（-10～10 m/s）分布圖、等效可燃?xì)庠企w積Q9隨時(shí)間變化圖。

圖7 不同風(fēng)速Q(mào)9隨時(shí)間變化值

由圖4可得風(fēng)速越大，天然氣云擴(kuò)散到達(dá)的最遠(yuǎn)距離越近，分別為55、45、40、37m。無(wú)風(fēng)狀態(tài)時(shí)，濃度沿泄漏源中軸線左右對(duì)稱分布，存在橫向風(fēng)力時(shí)，風(fēng)速越大，射流與泄漏源中軸線傾斜角度越大。射流核心區(qū)域濃度遠(yuǎn)高于邊界層，這是因?yàn)樘烊粴鈬娚涑龊?，其氣云邊界層不斷與大氣進(jìn)行質(zhì)量、動(dòng)量、能量交換，導(dǎo)致邊界層不斷卷吸空氣，濃度變低，風(fēng)速越大，邊界層卷吸空氣程度越強(qiáng)。同時(shí)，由于風(fēng)力的橫向輸送作用，風(fēng)速越大，氣云團(tuán)橫向輸送距離越遠(yuǎn)，射流中軸線傾斜角度越大?？偟膩?lái)看，4種風(fēng)速條件下射流中軸線傾斜角度均較小，這是因?yàn)槟M泄漏速率較大，動(dòng)量主導(dǎo)氣云團(tuán)運(yùn)動(dòng)軌跡［9］。

圖4 不同風(fēng)速XY切面氣云濃度

由圖5可得，速度矢量（VVEC）在V方向的分量分布圖與濃度分布圖類似，其邊界層值較小，射流中心區(qū)域值較大，這是因?yàn)闅庠茍F(tuán)邊界層不斷與大氣進(jìn)行動(dòng)量交換，導(dǎo)致速度衰減。又因?yàn)轱L(fēng)力的橫向輸送作用，速度云圖與泄漏源中軸線形成一定角度，風(fēng)速越大，角度越大。此外，同種風(fēng)速情況下泄漏源遠(yuǎn)端速度云圖傾斜角度大于近端，這是因?yàn)樘烊粴鈩傠x開泄漏孔時(shí)速度遠(yuǎn)大于風(fēng)速，風(fēng)速對(duì)其噴射影響很小，隨著徑向距離增大，天然氣速度逐漸衰減，大氣湍流的作用增強(qiáng)，天然氣速度等值線向右傾斜的趨勢(shì)越來(lái)越明顯［10］。

圖5 不同風(fēng)速XY切面氣云速度矢量在V方向的分量分布

由圖6可得，速度矢量（VVEC）在U方向的分量分布圖與濃度云圖不同，無(wú)風(fēng)時(shí)呈現(xiàn)“蝴蝶”形狀左右對(duì)稱分布，這是因?yàn)樯淞骱诵膮^(qū)域卷吸邊界層，導(dǎo)致對(duì)稱軸右側(cè)的天然氣產(chǎn)生向-X方向的速度，對(duì)稱軸左側(cè)的天然氣產(chǎn)生+X方向的速度。風(fēng)速逐漸增大時(shí)，監(jiān)測(cè)區(qū)域平面的氣云速度矢量在U方向的分量分布逐漸不規(guī)則，反映出湍流程度增強(qiáng)，但總體在射流中軸線左側(cè)的氣云速度大于右側(cè)氣云速度，這也是由于射流核心區(qū)域卷吸邊界層，導(dǎo)致射流軸線右側(cè)的天然氣產(chǎn)生向-X方向的速度，對(duì)稱軸左側(cè)的天然氣產(chǎn)生+X方向的速度。

圖6 不同風(fēng)速XY切面氣云速度矢量在U方向的分量

FLACS中用等效可燃?xì)庠企w積Q9來(lái)表示阻塞率低，通風(fēng)良好的開放場(chǎng)景與真實(shí)氣云等效化學(xué)計(jì)量比的均勻氣云體積，m3。Q9的計(jì)算公式為［9］：

式中，V為可燃?xì)怏w體積；BV為層流燃燒速度；E為在空氣中恒壓燃燒所產(chǎn)生的體積膨脹。

由圖7可得不同風(fēng)速條件下等效可燃?xì)庠企w積Q9隨時(shí)間變化趨勢(shì)基本一致。隨著泄漏開始，Q9在極短時(shí)間內(nèi)從0迅速增大，隨后增速放緩，最終達(dá)到最大值，此時(shí)泄漏仍在持續(xù)，但Q9只發(fā)生小幅度波動(dòng)，基本保持最大值不變。60 s時(shí)，隨著泄漏停止，Q9值迅速下降，在2 s之內(nèi)從最大值降為0。從縱向相比可以看出，在相同時(shí)刻，風(fēng)速越大，監(jiān)測(cè)區(qū)域內(nèi)的等效可燃?xì)庠企w積越小，這是因?yàn)殡S著風(fēng)速的增大，監(jiān)測(cè)區(qū)域內(nèi)氣體湍流程度增強(qiáng)，天然氣卷吸空氣程度增強(qiáng)，其濃度稀釋程度增強(qiáng)，導(dǎo)致Q9變小。此外，4種風(fēng)速條件下Q9值趨于穩(wěn)定后，風(fēng)速越大，Q9值波動(dòng)越大，這也是由于風(fēng)速越大，氣體湍流程度越強(qiáng)的原因造成的。

3.2 泄漏速率對(duì)氣云擴(kuò)散的影響

圖8—圖11分別為case1、case5、case6條件下泄漏后的天然氣擴(kuò)散達(dá)到穩(wěn)定狀態(tài)時(shí)Z=5 m的XY切面（泄漏源所在平面）天然氣體積分?jǐn)?shù)（1.5%～10%）分布圖、天然氣速度矢量（VVEC）在V方向上的分量（-35～35 m/s）分布圖、天然氣速度矢量在U方向上的分量（-10～10 m/s）分布圖、等效可燃?xì)庠企w積Q9隨時(shí)間變化圖。

由圖8可得，氣云圖沿泄漏源中軸線左右對(duì)稱分布，隨著泄漏速率增大，天然氣云團(tuán)擴(kuò)散到達(dá)的最遠(yuǎn)距離逐漸增大，分別為37、55、95 m。這是因?yàn)楫?dāng)泄漏孔徑一定（40 mm）時(shí)，泄漏速率越大，天然氣剛噴射出時(shí)初始動(dòng)能越大，在與大氣進(jìn)行質(zhì)量、動(dòng)量、能量交換過(guò)程中，其動(dòng)能衰減到0所需時(shí)間越長(zhǎng)，導(dǎo)致最終擴(kuò)散的最遠(yuǎn)距離越遠(yuǎn)。而且，泄漏速率越大，相同時(shí)間內(nèi)所噴射出的氣云體積越大，導(dǎo)致氣云覆蓋范圍越大。

圖8 不同泄漏速率XY切面氣云濃度

由圖9可得，與濃度分布云圖類似，速度矢量（VVEC）在V方向的分量分布（-35～35 m/s）云圖沿泄漏源中軸線左右對(duì)稱分布。隨著泄漏速率增大，V方向上的速度分量到達(dá)最遠(yuǎn)距離逐漸變大，分別為68、89、101 m，這是因?yàn)樾孤┛讖揭欢ǎ?0 mm）時(shí)，泄漏速率越大，剛噴射出的天然氣動(dòng)量越大，在氣云卷吸大氣的過(guò)程中，其速度衰減到0的過(guò)程越慢，導(dǎo)致速度邊界層到達(dá)距離越遠(yuǎn)。

圖9 不同泄漏速率XY切面氣云速度矢量在V方向的分量

由圖10得，速度矢量（VVEC）在U方向的分量（-10～10 m/s）分布圖均呈“蝴蝶”形狀左右對(duì)稱分布，這是因?yàn)樾孤┛滋幮纬韶?fù)壓區(qū)，射流不斷與大氣摻混造成的。隨著泄漏速率的增大，速度云圖分布面積增大。這是因?yàn)樾孤┧俾试酱?，射流核心區(qū)域?qū)Υ髿獾木砦饔迷綇?qiáng)造成的。除此之外，計(jì)算域邊界出現(xiàn)面積較大的速度云圖，這是由于射流噴射到最遠(yuǎn)端時(shí)，沿核心區(qū)域?qū)ΨQ軸形成“渦團(tuán)”，繼續(xù)卷吸大氣造成的，泄漏速率越大，“渦團(tuán)”尺寸越大［11］。

圖10 不同泄漏速率XY切面氣云速度矢量在U方向的分量分布

由圖11可得3種泄漏速率Q9變化趨勢(shì)基本相同，但從縱向相比可以看出，在相同時(shí)刻，泄漏速率越大，監(jiān)測(cè)區(qū)域內(nèi)的等效可燃?xì)庠企w積越大，這是因?yàn)樾孤┧俾试酱?，單位時(shí)間內(nèi)擴(kuò)散到濃度監(jiān)測(cè)區(qū)域的天然氣體積越大，導(dǎo)致Q9值越大。

圖11 不同泄漏速率Q9隨時(shí)間變化值

4 結(jié)論

1）通過(guò)對(duì)天然氣泄漏擴(kuò)散濃度分布和速度分布實(shí)驗(yàn)結(jié)果的數(shù)值模擬驗(yàn)證，證明FLACS數(shù)值模擬在研究天然氣泄漏擴(kuò)散濃度及速度分布規(guī)律具有一定可靠性；

2）橫向風(fēng)速越大，天然氣擴(kuò)散最遠(yuǎn)距離越近，濃度云圖向下風(fēng)向傾斜的趨勢(shì)越明顯；速度矢量在V方向的分量分布與濃度分布類似，風(fēng)速增大會(huì)使速度等值線向右傾斜的趨勢(shì)明顯；隨著橫向風(fēng)速增大，速度矢量在U方向的分量分布越不規(guī)則，但總體呈沿射流中軸線左側(cè)的氣云速度大于右側(cè)氣云速度；等效可燃?xì)庠企w積Q9趨于穩(wěn)定后值隨橫向風(fēng)速增大而減小，且波動(dòng)幅度增大；

3）泄漏速率越大，天然氣剛噴射出時(shí)初始動(dòng)能越大，速度衰減到0所需時(shí)間越長(zhǎng)，擴(kuò)散距離越遠(yuǎn)，穩(wěn)定后形成氣云覆蓋范圍越大，速度矢量在V方向的分量邊界層到達(dá)距離也越遠(yuǎn)；由于射流核心區(qū)域的卷吸大氣作用，射流中軸線左側(cè)形成天然氣產(chǎn)生+X方向速度，右側(cè)產(chǎn)生-X方向速度，速度矢量在U方向的分量分布云圖在泄漏源附近呈“蝴蝶”狀分布，泄漏速率越大，分布面積越大；射流末端出現(xiàn)“渦團(tuán)”，泄漏速率越大，“渦團(tuán)”尺寸越大；在相同時(shí)刻，泄漏速率越大，監(jiān)測(cè)區(qū)域內(nèi)的等效可燃?xì)庠企w積Q9越大。