基于FLACS的天然氣泄漏擴散數(shù)值模擬研究*

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

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表1 天然氣泄漏擴散工況

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

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

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

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

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

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

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

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

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

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

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

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

3.2 泄漏速率對氣云擴散的影響

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

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

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

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

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

由圖10得，速度矢量（VVEC）在U方向的分量（-10～10 m/s）分布圖均呈“蝴蝶”形狀左右對稱分布，這是因為泄漏孔處形成負壓區(qū)，射流不斷與大氣摻混造成的。隨著泄漏速率的增大，速度云圖分布面積增大。這是因為泄漏速率越大，射流核心區(qū)域?qū)Υ髿獾木砦饔迷綇娫斐傻?。除此之外，計算域邊界出現(xiàn)面積較大的速度云圖，這是由于射流噴射到最遠端時，沿核心區(qū)域?qū)ΨQ軸形成“渦團”，繼續(xù)卷吸大氣造成的，泄漏速率越大，“渦團”尺寸越大［11］。

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

由圖11可得3種泄漏速率Q9變化趨勢基本相同，但從縱向相比可以看出，在相同時刻，泄漏速率越大，監(jiān)測區(qū)域內(nèi)的等效可燃氣云體積越大，這是因為泄漏速率越大，單位時間內(nèi)擴散到濃度監(jiān)測區(qū)域的天然氣體積越大，導(dǎo)致Q9值越大。

圖11 不同泄漏速率Q9隨時間變化值

4 結(jié)論

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

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

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