燃氣泄漏過程及爆炸對建筑構件損毀的數(shù)值模擬研究

吳靜波，曾祥國*，楊 鑫，李 偉

(1.四川大學建筑與環(huán)境學院，深地科學與工程教育部重點實驗室，成都 610065；2.成都理工大學環(huán)境與土木工程學院，地質(zhì)災害防治與地質(zhì)環(huán)境保護國家重點實驗室，成都 610059)

近年來，由于燃氣泄漏引起爆炸事件屢屢發(fā)生，爆炸事故不僅會直接造成爆源附近的人員傷亡，而且會造成周圍建筑物損傷破壞，甚至會導致建筑物因某些關鍵結構構件的失效而喪失承載平衡能力，致使結構發(fā)生局部或整體的連續(xù)倒塌。因此，燃氣泄漏和擴散規(guī)律、燃氣泄漏引起爆炸的條件(如燃氣濃度、體積等)及其爆炸對建筑物損毀效應成為研究重點[1-4]。

目前，中外學者普遍單獨研究燃氣的泄漏過程或燃氣爆炸，針對燃氣泄漏后發(fā)生爆炸的兩個過程的綜合研究相對較少。現(xiàn)首先利用CFD軟件對管道燃氣泄漏到密閉空間和開放空間的規(guī)律進行模擬，接著根據(jù)泄漏時間和泄漏量對燃氣進行TNT當量等效換算，最后采用AUTODYN軟件對燃氣在密閉空間爆炸對三種不同類型的墻壁損毀效應進行了對比。以期為燃氣擴散后險情的控制、防護和救援提供依據(jù)。

1 管道燃氣泄漏過程CFD模擬

1.1 基本控制方程

在ANSYS/CFD模擬燃氣泄漏過程中，涉及的基本控制包括:連續(xù)性方程、動量方程、能量方程、組分方程和密度方程。

連續(xù)性方程:

(1)

式(1)中：ρ為流體密度；u為流體速度。

動量方程:

(2)

式(2)中：p為靜壓；τij為應力張量；gi和Fi分別為i方向上的重力體積力和外部體積力。

能量方程:

(3)

式(3)中:T表示流體溫度；k為流體的湍流導熱系數(shù)；σc為常數(shù)，一般取0.9～1；cpa為空氣的定壓比熱；cpv為泄漏物質(zhì)的定壓比熱；cp為混合物流體的定壓比熱，cp=ωcpv+(1-ω)cpa，ω為組分的質(zhì)量分率。

組分方程:當求解涉及化學反應或者組分問題時，CFD通過求解組分的對流擴散方程來預測每個組分的局部質(zhì)量分數(shù)，因研究的燃氣泄漏擴散涉及化學組分，所以還需要組分方程的控制，由組分質(zhì)量守恒定律可得氣體擴散組分方程[25]為

(4)

式(4)中：Dt為泄漏擴散流體的質(zhì)量擴散系數(shù)，Dt/D=(Sc/σc)μt/μ，σc是常數(shù)，一般情況為1；Sc為流體施密特數(shù)，Sc=μ/(ρD)。

密度方程:

(5)

式(5)中：R為普適氣體常數(shù)；P是絕對壓力；M為混合氣體的平均分子量；Ma為空氣的分子量；Mv為泄漏擴散物質(zhì)的分子量。

1.2 天然氣泄漏量的計算

管道泄漏模型分為小孔泄漏模型(泄漏點直徑小于20 mm)、管道泄漏模型(管道完全斷裂，泄漏孔面積為管道橫截面積的80%～100%)和其他泄漏模型[26]。管道因為施工開挖而開裂時，泄露口通常介于小孔和完全斷裂之間，綜合小孔泄漏模型和管道模型來計算泄漏量較為準確，根據(jù)文獻[27]知

(6)

式(6)中：一般對于燃氣速度較小、管道較長的情況，取n=1，而對于燃氣流速較大或者管道較短的情況，n=κ(κ為燃氣的等熵指數(shù))。如圖1所示，p1為管道起始端壓力，p2是泄漏點上游附近某點壓力，qm,0為泄漏點上方管道內(nèi)天然氣質(zhì)量流量，D是輸氣管道內(nèi)徑，λ是阻力系數(shù)，近似按城鎮(zhèn)燃氣設計規(guī)范提供的公式計算，L是泄漏點至管道起點的距離，取1 000 m。

圖1 天然氣管道泄漏示意圖

(7)

(8)

同時，作出假設

qm,0=qm

(9)

根據(jù)式(6)～式(9)，即可計算出泄漏量qm和泄漏處管內(nèi)壓力P3。表1為不同管徑和泄漏孔徑下天然氣泄漏量計算結果。

表1 天然氣泄漏量計算結果

1.3 計算模型

利用CFD 軟件對燃氣泄漏和擴散進行模擬分析，研究管道泄漏口附近泄漏擴散過程，分析計算域內(nèi)天然氣濃度變化規(guī)律和密閉空間壓力隨時間的變化規(guī)律。模擬考慮兩種工況，工況1為燃氣泄漏到密閉空間，工況2為燃氣泄漏到開放空間，其示意圖如圖2所示。依據(jù)中國現(xiàn)有的城鎮(zhèn)燃氣設計標準GB 50028—2006，并結合現(xiàn)有的城鎮(zhèn)燃氣管網(wǎng)鋪設概況，設燃氣管道運行壓力采用0.8 MPa，參考壓力為標準大氣壓，管道直徑為0.4 m，泄漏孔徑為0.2 m，其泄漏的質(zhì)量流量qm=15.875 kg/s。

圖2 兩種工況的簡化計算模型

1.4 數(shù)值模擬結果與分析

甲烷的爆炸極限是5%～15%(體積分數(shù))，因此主要分析天然氣泄漏擴散過程中其濃度場的變化。

如圖3所示，表示工況1不同時刻甲烷泄漏擴散的摩爾分數(shù)等值線圖。在泄漏口處天然氣的濃度最大，約為0.6。由于沒有外界因素的影響，天然氣泄漏擴散沿射流中心呈蘑菇狀對稱分布，在射流的中心線上天然氣的濃度最大，兩端逐漸降低。t=1 s時，從泄漏口向上至y=10 m，天然氣主要向上擴散，y=10～18 m之間，由于天然氣泄漏的初始速度的影響減小，蘑菇云迅速長大。t=2 s時，天然氣向上擴散到26 m，主要由于天然氣的密度比空氣小，其在浮力作用下向上運動。當t≥8 s時，天然氣擴散高度已經(jīng)達到計算域的上限，由于計算模型壁面的限制，擴散方向會發(fā)生改變，天然氣會沿著計算域的邊界向中心擴散，直至充滿整個計算域。

圖3 工況1不同時刻甲烷摩爾分數(shù)分布

圖4表示工況2不同時刻摩爾分數(shù)等值線，在t≤4 s時，天然氣的擴散形式和圖4結果類似，但其擴散速度較快。當擴散高度超過計算域時，計算域中只能顯示底部射流部分。

圖5表示監(jiān)測點位置，其中監(jiān)測點1位于泄漏口正上方1 m，對應圖6(a)黑色曲線。在擴散初期，甲烷的濃度先是急劇增大，達到0.63左右，然后快速下降，接著小幅度上升后趨于水平，同時可以看出其他四個點的濃度隨時間變化的趨勢與監(jiān)測點1類似；監(jiān)測點2和5處于同一高度，它們濃度變化起始時間十分接近，但是監(jiān)測點2濃度曲線增加到0.4后就基本趨于水平，而監(jiān)測點5濃度增加到0.27后卻迅速下降到0.1，其原因為監(jiān)測點2位于泄漏口的正上方，天然氣射流中心線附近流場比較穩(wěn)定，天然氣擴散到5 m(即監(jiān)測點2的高度)時，該監(jiān)測點濃度快速上升然后趨于不變；泄漏擴散剛發(fā)生不久，監(jiān)測點5所處位置的流場變化較大，以致其濃度變化范圍很大，當該區(qū)域流場趨于穩(wěn)定后，其變化就不再明顯。

圖5 監(jiān)測點位置

圖6 工況1與工況2監(jiān)測點甲烷質(zhì)量分數(shù)隨時間的變化曲線

結合圖6，可以發(fā)現(xiàn)兩種工況在0～15 s內(nèi)甲烷濃度曲線變化趨勢是類似的，隨后工況1甲烷濃度由于空間的限制會緩慢地上升，而工況2因為擴散到開放空間，天然氣能夠自由流動，所以監(jiān)測點的濃度曲線變化不明顯。

上述分析僅通過某些監(jiān)測點的質(zhì)量分數(shù)來描述天然氣泄漏的狀態(tài)，但不足以反映整個流場甲烷濃度隨時間的變化情況。圖7為工況1甲烷質(zhì)量和摩爾分數(shù)隨時間的變化曲線，可知在10 s前，兩條曲線的增長速率是比較緩慢的，隨后以較大的斜率增加；甲烷摩爾分數(shù)比質(zhì)量分數(shù)大，且兩條曲線基本上呈線性增長，這與實際情況相符。圖8為工況1計算域內(nèi)平均壓力(相對標準大氣壓而言)隨甲烷泄漏時間的變化規(guī)律，隨著甲烷擴散的進行，其先快速地增加，并在60 s后逐漸穩(wěn)定。通過擬合，可知平均壓力隨時間呈現(xiàn)3次多項式變化趨勢。

圖7 工況1甲烷質(zhì)量和摩爾分數(shù)隨時間的變化曲線

圖7和圖8中的曲線分別充分地反映了泄漏后流場內(nèi)甲烷濃度和壓力的變化，這對實時了解天然氣的濃度分布和泄漏域內(nèi)壓力大小，以及泄漏擴散后險情的控制、防護和救援具有重要科學意義。

圖8 工況1空間內(nèi)平均壓力隨時間的變化曲線

2 燃氣爆炸對墻壁損毀的數(shù)值模擬

2.1 TNT當量計算

關于可燃氣體與TNT當量法之間換算與適用性可參考文獻[14]。

甲烷/空氣混合氣體的爆炸化學反應方程式為

(10)

氧氣在空氣中的體積濃度為21%，1 mol甲烷在空氣中完全反應的化學方程式為

(11)

根據(jù)蓋斯定律，參考表2甲烷/空氣混合氣體反應生成焓，可以求出1 mol乙炔在空氣中反應的爆熱為Qv，則單位質(zhì)量的甲烷/空氣混合氣體爆炸產(chǎn)生的爆熱Qm為

表2 甲烷/空氣混合氣體及爆炸產(chǎn)物生成焓

(12)

甲烷泄漏質(zhì)量為

Mm=qmt

(13)

式(13)中：Mm為甲烷泄漏質(zhì)量；qm為甲烷泄漏的質(zhì)量流量，取15.875 kg/s；t為泄漏時間，21.7 s。

等效TNT質(zhì)量為

(14)

式(14)中：MTNT為等效TNT質(zhì)量；α為TNT當量效率，取0.04；Qm為甲烷/空氣混合氣體的爆熱；QTNT為TNT的爆熱，取4.52 MJ/kg。

根據(jù)式(12)～式(14)，計算獲得:Qm=3 283.55 KJ/kg，Mm=344.49 kg，MTNT=10 kg。

2.2 爆炸對墻壁損毀的數(shù)值模擬結果與探討

選擇1.3工況1(密閉空間)作為模擬背景，以磚墻、鋼筋混凝土墻和混凝土墻作為密閉空間不同種類的墻壁。數(shù)值模擬采用Autodyn軟件，計算中涉及的材料包括:空氣、TNT、混凝土、鋼筋、磚和砂漿，各材料模型見表3。

表3 材料模型

考慮到計算量的大小，數(shù)值模擬過程中選擇4 m×1.8 m×0.1 m作為墻壁計算模型，圖9為1/2模型。鋼筋與混凝土采用分離式建模，以共節(jié)點方式耦合，邊緣鋼筋距混凝土模型上下距離為0.2 m、右側為0.4 m，鋼筋網(wǎng)間距為0.2 m×0.2 m，橫截面積7.85×10-4m2；磚塊尺寸為0.225 m×0.075 m×0.1 m。固定墻壁最下端，覆蓋墻壁的空氣采用Flow-out邊界條件；計算時間為10 ms。以1/2模型左下角為原點，點1坐標(x，y，z)=(0.4，0.2，0)和點7坐標(x，y，z)=(0，0.6，0)，相鄰點間距0.2 m，所有觀察點均位于墻壁的迎爆面上且三種墻壁的觀察點分布一致。此外，參考根據(jù)1.4節(jié)工況1模擬結果，選取爆源中心到墻壁距離為0.5 m。為避免計算時間過多地消耗在流體單元，實現(xiàn)遠場爆炸，基于Autodyn的remap技術，將一維計算結果映射到三維計算模擬中，其映射部位在左下角正前方。

圖9 磚墻和鋼筋混凝土墻模型及觀察點分布

圖10表示三種類型墻壁爆炸后的效果。可以看出，三種類型墻壁均在其下端中心部位產(chǎn)生凹陷變形或破壞，這是因爆炸沖擊作用向外擴展造成的。圖10(a)中磚墻中下部因爆炸強沖擊作用而產(chǎn)生方形斷口，其面積約為0.076 m2，最下方磚塊也產(chǎn)生一條破壞帶，破壞長度為1.2 m。圖10(b)和圖10(c)分別表示鋼筋混凝土墻和混凝土墻的損傷情況，兩墻壁破壞部位基本相似，墻壁上側和左右側產(chǎn)生多條層裂損傷，這是反射拉伸應力波在自由面附近造成的。同時兩墻壁都在底部位產(chǎn)生兩條破壞帶，前者最底部位和最底偏上部位破壞帶長度分別約為1.09 m和0.56 m，后者最底部位和最底偏上部位破壞帶長度分別約為1.05 m和0.2 m。

圖10 三種類型墻爆后效果

如圖11(a)～圖11(c)，表示三種墻壁壓力在1 ms內(nèi)的曲線。可以看出，磚墻受到的爆炸壓力衰減速度快于鋼筋混凝土墻和混凝土墻所受到的壓力，后兩者壓力出現(xiàn)明顯的負壓，且后期壓力震蕩性較大。圖11(d)中以觀察點Gauge#1和#7分別作為水平和豎直方向上觀察點的x軸坐標起點(下同)，表示各觀察點壓力峰值隨距離變化的曲線。除磚墻觀察點Gauge#7～12壓力峰值隨距離變化表現(xiàn)波動外，其余各觀察點壓力峰值隨距離增大而下降。爆炸應力波正向接觸墻壁后，墻內(nèi)壓縮應力波便向四周傳播，并隨傳播距離而發(fā)生能量衰減，由于應力波在墻壁的豎向先于水平發(fā)生反射，并與壓縮應力波發(fā)生疊加，從而降低墻壁豎向壓力峰值。故三種墻壁水平觀察點Gauge#1～6的壓力峰值多大于豎向觀察點Gauge#7～12的壓力峰值。該現(xiàn)象在距離為0～0.6 m范圍內(nèi)表現(xiàn)明顯，而距離為0.6～1 m范圍橫豎向壓力峰值較接近。圖11(e)為壓力峰值隨距離變化的擬合曲線。除磚墻觀察點Gauge#7～12之外，其余觀察點都表現(xiàn)為指數(shù)衰減趨勢，且水平觀察點壓力峰值衰減幅度大于豎向觀察點。在相同距離下，磚墻橫豎向壓力峰值的差值大于其余兩種墻壁，鋼筋混凝土墻橫豎向壓力峰值的差值最小。

圖11 壓力時程曲線與壓力峰值和距離的關系

如圖12(a)～圖12(c)，表示三種墻壁速度時程曲線。由于受到壓縮應力波沖擊作用，墻壁壓力陡然上升，造成墻壁質(zhì)點速度也迅速增大到峰值，可知0～1 ms內(nèi)加速度值很大；隨著時間增加，墻壁受力必然下降，而后期的速度達到峰值后卻是緩慢下降，幾乎是呈直線下降趨勢，表明加速度變化很小。此外，還可看出磚墻的觀察點速度大于鋼筋混凝土墻和混凝土墻上對應的觀察點速度。在圖12(d)中，隨著距離增大，三種墻壁觀察點的速度峰值都快速下降；相同距離下磚墻的觀察點速度峰值都大于鋼筋混凝土墻和混凝土墻的觀察點速度峰值，而后兩者的速度峰值較接近。相比鋼筋混凝土墻與混凝土墻，磚和砂漿的動態(tài)抗壓強度都較小，且砂漿更容易受到破壞而喪失抵抗能力，造成磚墻抗爆炸壓力較弱，墻壁速度較大而更易破壞。圖12(e)為磚墻和鋼筋混凝土墻的觀察點速度峰值的擬合曲線。由于鋼筋混凝土墻和混凝土墻的觀察點速度峰值接近且變化趨勢一致，因此采用鋼筋混凝土墻觀察點速度峰值擬合曲線代替后者的擬合曲線。通過對磚墻和鋼筋混凝土墻的水平與豎直觀察點的速度峰值擬合，發(fā)現(xiàn)這兩類墻壁觀察點速度峰呈現(xiàn)指數(shù)衰減趨勢，且磚墻觀察點速度峰值衰減幅度更大。

圖12 速度時程曲線與速度峰值和距離的關系

如圖13(a)～圖13(c)表示三種墻壁測點位移時程曲線。可以看出，越靠近爆源，位移越大；對磚墻而言，Gauge#7的位移從墻壁產(chǎn)生位移開始就大于Gauge#1的位移，而鋼筋混凝土墻和混凝土墻則大約在2.5 ms后Gauge#7的位移才大于Gauge#1的位移。圖13(d)為磚墻和鋼筋混凝土墻Gauge#7的位移時程曲線，從爆炸應力波加載墻壁開始，磚墻的Gauge#7的位移就明顯大于鋼筋混凝土墻的Gauge#7的位移，這是因磚墻抗壓強度較小，更易于發(fā)生變形和破壞。通過擬合，發(fā)現(xiàn)兩種墻壁相同位置觀察點都呈線性關系。

圖13 位移時程曲線與位移峰值和距離的關系

圖13(e)為位移峰值與距離的關系。三種墻壁豎直觀察點Gauge#7～12的位移峰值基本上都大于水平觀察點Gauge#1～6，且相同距離下磚墻的豎直觀察點與水平觀察點位移峰值的差值大于其余兩種墻壁的差值。由于墻壁下端固定，爆炸應力波正面作用墻壁，使墻壁朝后發(fā)生變形，而墻壁兩側因下端固定而變形相對較小，因此墻壁豎直觀測點位移大于水平觀測點位移。圖13(f)表示位移峰值擬合曲線，這兩類墻壁觀察點位移峰同樣呈現(xiàn)指數(shù)衰減趨勢。

3 結論

(1)對比燃氣泄漏到密閉空間與開放空間兩種工況的濃度分布圖，發(fā)現(xiàn)前者擴散速率慢于后者，主要由于密閉空間的限制降低了甲烷和空氣的對流擴散速率。

(2)兩種工況在擴散時間為0～15 s內(nèi)，甲烷質(zhì)量分數(shù)變化曲線基本一致；擴散時間超過15 s后，密閉空間的甲烷質(zhì)量分數(shù)會繼續(xù)增大，而開放空間的甲烷質(zhì)量分數(shù)則保持基本不變，這表明甲烷泄漏擴散時間對其質(zhì)量分數(shù)有重要影響。

(3)三種墻壁觀察點壓力峰值、速度峰值、位移峰值基本上隨距離增大而表現(xiàn)為指數(shù)衰減趨勢。磚墻水平壓力峰值基本上大于豎直壓力峰值，水平速度峰值、位移峰值卻分別小于豎直速度峰值、位移峰值；鋼筋混凝土墻和混凝土墻皆表現(xiàn)為水平壓力峰值、位移峰值小于豎直壓力峰值、位移峰值，水平速度峰值大于豎直速度峰值。

(4)通過CFD計算獲得了甲烷平均濃度和壓力隨時間變化規(guī)律，同時運用Autodyn模擬燃氣爆炸對密閉空間墻壁損毀效應，兩者共同為甲烷泄漏擴散后險情的控制、防護、救援以及對建筑物損毀的檢測、評估提供了科學依據(jù)。