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

(1.四川大學(xué)建筑與環(huán)境學(xué)院，深地科學(xué)與工程教育部重點(diǎn)實(shí)驗(yàn)室，成都 610065；2.成都理工大學(xué)環(huán)境與土木工程學(xué)院，地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點(diǎn)實(shí)驗(yàn)室，成都 610059)

近年來，由于燃?xì)庑孤┮鸨ㄊ录覍野l(fā)生，爆炸事故不僅會(huì)直接造成爆源附近的人員傷亡，而且會(huì)造成周圍建筑物損傷破壞，甚至?xí)?dǎo)致建筑物因某些關(guān)鍵結(jié)構(gòu)構(gòu)件的失效而喪失承載平衡能力，致使結(jié)構(gòu)發(fā)生局部或整體的連續(xù)倒塌。因此，燃?xì)庑孤┖蛿U(kuò)散規(guī)律、燃?xì)庑孤┮鸨ǖ臈l件(如燃?xì)鉂舛取Ⅲw積等)及其爆炸對建筑物損毀效應(yīng)成為研究重點(diǎn)[1-4]。

目前，中外學(xué)者普遍單獨(dú)研究燃?xì)獾男孤┻^程或燃?xì)獗ǎ槍θ細(xì)庑孤┖蟀l(fā)生爆炸的兩個(gè)過程的綜合研究相對較少?，F(xiàn)首先利用CFD軟件對管道燃?xì)庑孤┑矫荛]空間和開放空間的規(guī)律進(jìn)行模擬，接著根據(jù)泄漏時(shí)間和泄漏量對燃?xì)膺M(jìn)行TNT當(dāng)量等效換算，最后采用AUTODYN軟件對燃?xì)庠诿荛]空間爆炸對三種不同類型的墻壁損毀效應(yīng)進(jìn)行了對比。以期為燃?xì)鈹U(kuò)散后險(xiǎn)情的控制、防護(hù)和救援提供依據(jù)。

1 管道燃?xì)庑孤┻^程CFD模擬

1.1 基本控制方程

在ANSYS/CFD模擬燃?xì)庑孤┻^程中，涉及的基本控制包括:連續(xù)性方程、動(dòng)量方程、能量方程、組分方程和密度方程。

連續(xù)性方程:

(1)

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

動(dòng)量方程:

(2)

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

能量方程:

(3)

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

組分方程:當(dāng)求解涉及化學(xué)反應(yīng)或者組分問題時(shí)，CFD通過求解組分的對流擴(kuò)散方程來預(yù)測每個(gè)組分的局部質(zhì)量分?jǐn)?shù)，因研究的燃?xì)庑孤U(kuò)散涉及化學(xué)組分，所以還需要組分方程的控制，由組分質(zhì)量守恒定律可得氣體擴(kuò)散組分方程[25]為

(4)

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

密度方程:

(5)

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

1.2 天然氣泄漏量的計(jì)算

管道泄漏模型分為小孔泄漏模型(泄漏點(diǎn)直徑小于20 mm)、管道泄漏模型(管道完全斷裂，泄漏孔面積為管道橫截面積的80%～100%)和其他泄漏模型[26]。管道因?yàn)槭┕ら_挖而開裂時(shí)，泄露口通常介于小孔和完全斷裂之間，綜合小孔泄漏模型和管道模型來計(jì)算泄漏量較為準(zhǔn)確，根據(jù)文獻(xiàn)[27]知

(6)

式(6)中：一般對于燃?xì)馑俣容^小、管道較長的情況，取n=1，而對于燃?xì)饬魉佥^大或者管道較短的情況，n=κ(κ為燃?xì)獾牡褥刂笖?shù))。如圖1所示，p1為管道起始端壓力，p2是泄漏點(diǎn)上游附近某點(diǎn)壓力，qm,0為泄漏點(diǎn)上方管道內(nèi)天然氣質(zhì)量流量，D是輸氣管道內(nèi)徑，λ是阻力系數(shù)，近似按城鎮(zhèn)燃?xì)庠O(shè)計(jì)規(guī)范提供的公式計(jì)算，L是泄漏點(diǎn)至管道起點(diǎn)的距離，取1 000 m。

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

(7)

(8)

同時(shí)，作出假設(shè)

qm,0=qm

(9)

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

表1 天然氣泄漏量計(jì)算結(jié)果

1.3 計(jì)算模型

利用CFD 軟件對燃?xì)庑孤┖蛿U(kuò)散進(jìn)行模擬分析，研究管道泄漏口附近泄漏擴(kuò)散過程，分析計(jì)算域內(nèi)天然氣濃度變化規(guī)律和密閉空間壓力隨時(shí)間的變化規(guī)律。模擬考慮兩種工況，工況1為燃?xì)庑孤┑矫荛]空間，工況2為燃?xì)庑孤┑介_放空間，其示意圖如圖2所示。依據(jù)中國現(xiàn)有的城鎮(zhèn)燃?xì)庠O(shè)計(jì)標(biāo)準(zhǔn)GB 50028—2006，并結(jié)合現(xiàn)有的城鎮(zhèn)燃?xì)夤芫W(wǎng)鋪設(shè)概況，設(shè)燃?xì)夤艿肋\(yùn)行壓力采用0.8 MPa，參考壓力為標(biāo)準(zhǔn)大氣壓，管道直徑為0.4 m，泄漏孔徑為0.2 m，其泄漏的質(zhì)量流量qm=15.875 kg/s。

圖2 兩種工況的簡化計(jì)算模型

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

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

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

圖3 工況1不同時(shí)刻甲烷摩爾分?jǐn)?shù)分布

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

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

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

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

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

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

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

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

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

2 燃?xì)獗▽Ρ趽p毀的數(shù)值模擬

2.1 TNT當(dāng)量計(jì)算

關(guān)于可燃?xì)怏w與TNT當(dāng)量法之間換算與適用性可參考文獻(xiàn)[14]。

甲烷/空氣混合氣體的爆炸化學(xué)反應(yīng)方程式為

(10)

氧氣在空氣中的體積濃度為21%，1 mol甲烷在空氣中完全反應(yīng)的化學(xué)方程式為

(11)

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

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

(12)

甲烷泄漏質(zhì)量為

Mm=qmt

(13)

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

等效TNT質(zhì)量為

(14)

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

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

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

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

表3 材料模型

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

圖9 磚墻和鋼筋混凝土墻模型及觀察點(diǎn)分布

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

圖10 三種類型墻爆后效果

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

圖11 壓力時(shí)程曲線與壓力峰值和距離的關(guān)系

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

圖12 速度時(shí)程曲線與速度峰值和距離的關(guān)系

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

圖13 位移時(shí)程曲線與位移峰值和距離的關(guān)系

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

3 結(jié)論

(1)對比燃?xì)庑孤┑矫荛]空間與開放空間兩種工況的濃度分布圖，發(fā)現(xiàn)前者擴(kuò)散速率慢于后者，主要由于密閉空間的限制降低了甲烷和空氣的對流擴(kuò)散速率。

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

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

(4)通過CFD計(jì)算獲得了甲烷平均濃度和壓力隨時(shí)間變化規(guī)律，同時(shí)運(yùn)用Autodyn模擬燃?xì)獗▽γ荛]空間墻壁損毀效應(yīng)，兩者共同為甲烷泄漏擴(kuò)散后險(xiǎn)情的控制、防護(hù)、救援以及對建筑物損毀的檢測、評估提供了科學(xué)依據(jù)。