王桂林，歐陽嘯天，翟 俊，孫 帆

（1. 重慶大學土木工程學院，重慶 400045；2. 庫區(qū)環(huán)境地質(zhì)災害防治國家地方聯(lián)合工程研究中心，重慶 400045；3. 重慶大學三峽庫區(qū)生態(tài)環(huán)境教育部重點實驗室，重慶 400045）

城市地下綜合管廊含有燃氣管道，管道燃氣泄漏爆炸事故時有發(fā)生，如2014 年中國臺灣高雄燃氣管廊爆炸[1]（圖1），2016 年德國綜合管廊燃氣爆炸等[2]。由于管廊常位于人員密集區(qū)，地下管廊一旦發(fā)生爆炸，將會給地面人員的生命和財產(chǎn)造成巨大損失，因此，研究地下管廊燃氣爆炸地面響應對安全生產(chǎn)和生活具有重要意義。

由于實驗成本高且安全風險大，針對地下空間爆炸問題，學者們一般采用有限元數(shù)值模擬方法[3]開展研究。曲樹盛等[4]利用AUTODYN 有限元軟件模擬了地鐵車站內(nèi)爆炸波的傳播過程和衰減規(guī)律，廖維張等[5-6]利用ANSYS/LS-DYNA 有限元軟件模擬了爆炸波在二層地鐵站的傳播規(guī)律和中柱的動力響應，田力等[7]利用ABAQUS 和LS-DYNA 軟件分析了地下空間爆炸產(chǎn)生的荷載對建筑物的影響。然而在解決由爆炸引起的巖土體大變形問題方面，有限元方法常采用輸入應力波的方式模擬爆炸應力，這與氣體爆炸沖擊巖土體的實際情況存在較大差異，此外，巖土體大變形將造成有限元網(wǎng)格畸變。物質(zhì)點法采用拉格朗日質(zhì)點和歐拉網(wǎng)格雙重描述，將連續(xù)體離散成一組質(zhì)點，每個質(zhì)點代表一塊材料區(qū)域并攜帶該材料區(qū)域的所有物理信息，計算網(wǎng)格僅用于動量方程的求解和空間導數(shù)的計算[8]，同時不需要對物質(zhì)交界面進行特殊處理，自動滿足無滑移的接觸條件[9]，能夠直接模擬物體之間的沖擊現(xiàn)象。這些特點使得物質(zhì)點法在處理爆炸沖擊[10-12]和巖土大變形問題上[13-14]較有限元方法更有優(yōu)勢。

圖1 2014 年高雄地下丙烯管道爆炸現(xiàn)場[1]Fig. 1 Explosion scene of the underground propylene pipeline in Kaohsiung in 2014[1]

本研究依托重慶市某地下綜合管廊試點工程，采用點火增長模型，運用物質(zhì)點法模擬淺埋管廊甲烷氣體泄漏后爆炸沖擊管廊本體結(jié)構(gòu)和圍巖的過程，以獲得地面壓強與位移的響應特性。

1 模型設(shè)置

1.1 接觸算法

標準格式的物質(zhì)點法可以自動處理物體之間的無滑移接觸[8]，但為考慮物體之間的分離和滑動，Bardenhagen 等[15]將接觸算法引入物質(zhì)點法中，通過在網(wǎng)格節(jié)點上引入多重速度場的方式，模擬物體之間由于速度不同而產(chǎn)生的分離和滑動。該算法中，接觸判據(jù)表達式為

如果物體B 與C 之間發(fā)生接觸和穿透，則需要施加接觸力以防止穿透發(fā)生，物體C 對B 產(chǎn)生的接觸力計算公式為

通過在單步計算中遍歷對單個網(wǎng)格節(jié)點上產(chǎn)生影響的多個物體間的兩兩接觸，將得到的接觸力進行矢量疊加，從而實現(xiàn)處理多個物體之間接觸問題。計算示意圖如圖2 所示。

圖2 多物體接觸計算示意圖Fig. 2 Schematic diagram of multi-object contact algorithm

1.2 甲烷爆轟模型設(shè)置

甲烷的燃燒有定壓燃燒、爆燃、爆轟等形式，其中爆轟最劇烈，且對周圍環(huán)境的破壞性最強。本研究采用點火增長模型模擬管廊中泄漏甲烷氣體的爆轟過程，并將兩種常見的定常爆轟反應率函數(shù)結(jié)合使用。

1.2.1 反應率函數(shù)設(shè)定

Wilkins 函數(shù)表達式為

式中：ps為反應完全時爆炸產(chǎn)物壓強函數(shù)，E0為氣體的單位體積內(nèi)能。

1.2.2 人工體積黏性

在爆轟過程中沖擊波前后的物質(zhì)點參數(shù)變化巨大，在應力項上施加人工體積黏性項q[16]，將爆轟波的強間斷面化成一個急劇但連續(xù)變化的波陣面，起到光滑間斷面的作用。

人工體積黏性計算公式表示為

式中：c0、c1為輸入常數(shù)，分別取 3.0、0.3；c 為聲速； ρ為密度；為體積應變率。

1.3 模型及參數(shù)

重慶市某地下綜合管廊試點工程的綜合管廊標準斷面高3.3 m，寬9.0 m，管廊包括燃氣艙、電力艙和水信艙3 個艙，如圖3 所示，頂板位于途經(jīng)道路設(shè)計路面標高以下約3.4 m。建立的模型如圖4 所示。

圖3 某地下綜合管廊現(xiàn)場照片F(xiàn)ig. 3 Picture of an underground comprehensive pipe gallery

圖4 數(shù)值模擬模型內(nèi)部示意圖Fig. 4 Internal schematic diagram of numerical simulation model

模型尺寸為40.0 m × 30.0 m × 17.0 m，質(zhì)點均勻分布，質(zhì)點間距0.5 m，正方體網(wǎng)格邊長1.0 m，質(zhì)點總數(shù)為157 632 個。燃氣艙、電力艙和水信艙的截面尺寸分別為2.0 m × 3.0 m、4.0 m × 3.0 m、2.0 m × 3.0 m，管廊結(jié)構(gòu)厚度均為0.5 m。由于實際情況中管廊上覆巖土體與管廊之間的黏結(jié)力小，故模擬中不考慮。參考Zhu 等[17]開展的混合氣體爆炸模擬實驗，設(shè)置發(fā)生泄漏的甲烷-空氣混合氣體充滿管廊中部，長度為4.0 m 的區(qū)域，起爆點位于混合氣體正中心。模型四周采用黏彈性邊界模擬側(cè)向無限巖土體，底部采用對稱邊界模擬基巖面，頂部采用自由邊界模擬地面。

模型圍巖和管廊本體結(jié)構(gòu)模擬采用DP 強度準則，甲烷-空氣混合氣體模擬采用高能燃燒模型。參考工程地質(zhì)勘察報告和文獻[18]，巖土物理力學參數(shù)取值如表1 所示，混合氣體參數(shù)取值如表2 所示。其中：E 為彈性模量， ν為泊松比， σt為 巖石抗拉強度， qφ、 Kφ為由摩爾庫倫參數(shù)換算得到的材料參數(shù)，ψ為剪脹角，Vf為體積分數(shù)。

表1 圍巖和管廊本體結(jié)構(gòu)力學參數(shù)Table 1 Mechanical parameters of surrounding rock and pipe gallery structure

表2 甲烷-空氣混合氣體高能燃燒模型計算參數(shù)Table 2 Calculation parameters of high-energy combustion model of methane-air mixture

2 模型驗證

參照 Zhu 等[17]、Xu 等[19]和 Qu 等[20]的實驗，在一個長900.0 m、一端封閉、截面面積為7.2 m2的正方形隧道中，充滿區(qū)域長度約為14 m、體積為100 m3、體積分數(shù)為10%的甲烷-空氣混合氣體。點火點在封閉隧道的一端。圖5 為物質(zhì)點法模型示意圖，其中質(zhì)點間距0.3 m，網(wǎng)格間距0.6 m，質(zhì)點總數(shù)363 001。表3 為數(shù)值模擬結(jié)果及其與文獻實驗結(jié)果的對比，其中： δ為相對誤差。

圖5 模型截面示意圖Fig. 5 Model section schematic

表3 數(shù)值模擬驗證結(jié)果Table 3 Numerical simulation verification results

分析表3 可以發(fā)現(xiàn)，物質(zhì)點法數(shù)值模擬結(jié)果與實驗結(jié)果的偏差約為9%，較能反映實際情況[21]，產(chǎn)生誤差的原因主要有以下3 點：（1）由于物質(zhì)點法的算法特性，單個測量質(zhì)點不能代表具體的物質(zhì)形狀，相較于實驗中測量元件具有明確的尺寸，物質(zhì)點法中測量質(zhì)點得出的結(jié)果與實際情況存在出入；（2）本研究中質(zhì)點間距為0.3 m，網(wǎng)格邊長為0.6 m，對于精確模擬來說精度不夠，然而隨著物質(zhì)點數(shù)量增加，網(wǎng)格劃分密度增大，模擬結(jié)果會更加貼近實際情況；（3）實驗中隧道與高速氣體存在一定的摩擦力，這種摩擦力會影響爆炸壓強，而具體區(qū)域的摩擦力大小受氣體的流速和隧道壁的粗糙程度影響，致使數(shù)值模擬結(jié)果大于實際情況。

3 爆炸作用地面響應特性

3.1 地面壓強響應特性

3.1.1 地面響應壓強平面分布特性

圖6 為混合氣體引爆0.02、0.20、0.50 s 時對應的地面壓強分布情況。由圖6 可以看出，起爆后0.02 s，地面上以起爆點豎直方向為中心的橢圓形區(qū)域出現(xiàn)方向向上的響應壓強，并向四周輻射，這是由于混合氣體引爆以后，受到爆轟壓力的管廊襯砌產(chǎn)生向上的形變，從而沖擊地面造成的。起爆后0.20 s，爆炸應力波向四周傳播，管廊上覆地面響應壓強方向整體向下，這是由于地表為自由邊界，應力波傳播至地表后發(fā)生反射導致向下傳播；部分區(qū)域出現(xiàn)應力集中現(xiàn)象，是由于應力波在管廊水信艙及電力艙空洞中發(fā)生了反射；緊鄰管廊上覆地表的巖土體出現(xiàn)應力集中現(xiàn)象，是由于應力波在通過管廊和巖土接觸截面時發(fā)生了折射。觀察地面上整體響應壓強分布情況，可以發(fā)現(xiàn)燃氣艙一側(cè)地面響應壓強分布集中，數(shù)值較大；電力艙一側(cè)響應壓強分布范圍更廣，但數(shù)值較小。起爆后0.50 s，地面響應壓強分布逐漸穩(wěn)定，在模型四周出現(xiàn)的響應壓強是由于吸收邊界未能完全吸收反射波。地面中部出現(xiàn)的應力集中現(xiàn)象是由于管廊未考慮阻尼，在管廊內(nèi)部反射的應力波未完全衰減。

圖6 爆炸后不同時刻地面壓強分布Fig. 6 Distribution of ground pressure at different times after explosion

3.1.2 管廊橫向壓強分布特性

圖7 為管廊橫向距起爆點水平距離-4.0、0、4.0、8.0 m 處質(zhì)點壓強隨時間的變化曲線。由圖7可知，從響應時間來看，地面質(zhì)點的響應時間隨距起爆點水平距離的增大而延長。管廊上覆地面質(zhì)點（0、4.0、8.0 m）在經(jīng)歷初始響應壓強峰值后均會出現(xiàn)方向向下的負壓強，而圍巖上覆地面（-4.0、-8.0 m）則不會出現(xiàn)這種情況。此外，除了起爆點正上方質(zhì)點外，其他質(zhì)點在爆炸過程中出現(xiàn)了由于接觸界面反射作用形成的次生響應壓強波峰，且隨著距離的增大，次生波峰峰值逐漸增大。根據(jù)主波和次生波產(chǎn)生的階段，可以將爆炸響應分成3 個階段：爆炸響應階段、次生響應階段和漸趨穩(wěn)定階段。

圖8 為管廊橫向質(zhì)點壓強峰值隨起爆點距離的變化曲線。由圖8 可以看出，爆炸引起管廊橫向的地面峰值響應壓強受地下情況影響明顯。其中，響應壓強峰值的最大值偏離起爆點正上方，管廊上覆地面距起爆點2.0 m 處峰值壓強最大，為37.16 kPa。從響應壓強峰值分布來看，圍巖上覆地面響應峰值的最大值小于管廊上覆地面，但隨距起爆點距離的增大，其減小幅度小于管廊上覆地面。響應壓強最大值峰值點的偏移，是由于淺埋管廊與巖土體接觸面和艙體臨空面的存在，使得爆炸應力波傳輸?shù)焦芾缺倔w結(jié)構(gòu)時發(fā)生折射，導致響應壓強峰值出現(xiàn)偏移現(xiàn)象。

圖7 管廊橫向質(zhì)點壓強隨時間變化曲線Fig. 7 Variation curves of particle pressure in the horizontal direction of pipe gallery with time

圖8 管廊橫向質(zhì)點壓強峰值隨起爆點距離的變化Fig. 8 Peak particle pressure in the horizontal direction of pipe gallery with distance of initiation point

3.1.3 管廊縱向壓強分布特性

由于模型在縱向?qū)ΨQ，取縱向距起爆點水平距離0、2.0、4.0、6.0 m 處質(zhì)點的響應壓強隨時間的變化曲線，如圖9 所示。由圖9 可知，地面質(zhì)點的響應時間隨距起爆點水平距離的增大而延長，但與橫向的響應特性不同，縱向各質(zhì)點在受爆炸作用后產(chǎn)生的次生波振幅較小，且隨距離的變化不大；此外，處于管廊正上方縱向各質(zhì)點的壓強均豎直向下。

圖10 為管廊縱向質(zhì)點壓強峰值隨起爆點距離的變化曲線。從圖10 中可以看出，沿管廊縱向的地面響應峰值壓強隨與起爆點水平距離的增大而經(jīng)歷先增大后減小的過程，且隨著距起爆點距離的增大，壓強峰值變化速率逐漸減小。最大值出現(xiàn)在距管廊起爆點2.0 m 處，峰值壓強為26.29 kPa。出現(xiàn)這種現(xiàn)象的原因是：混合氣體在橫向尺寸較大，氣體發(fā)生化學反應的先后不同導致應力波疊加，出現(xiàn)峰值壓強最大值偏離起爆點；距離起爆點越遠，混合氣體中甲烷消耗得越多，釋放的壓強越小，因而對距離起爆點較遠的地面質(zhì)點影響更小。

圖9 管廊縱向質(zhì)點壓強隨時間變化曲線Fig. 9 Variation curves of particle pressure in the longitudinal direction of the pipe gallery with time

圖10 管廊縱向質(zhì)點壓強峰值隨起爆點距離變化Fig. 10 Peak particle pressure of the pipe gallery with the distance of the initiation point in the longitudinal direction

3.2 地面豎向位移響應特性

3.2.1 地面響應豎向位移分布特性

圖11 為截取地表各質(zhì)點爆炸作用下的位移情況。（1）從絕對位移來看，模型地面整體出現(xiàn)了豎直向下的位移，這是由于管廊可燃性氣體爆炸時，應力波向四周圍巖傳播，管廊底部圍巖受力產(chǎn)生壓縮變形，使得地面發(fā)生沉降。（2）從地面隆起形態(tài)來看，在隆起范圍方面，由于管廊燃氣艙中甲烷爆炸作用，管廊燃氣艙上覆地面立即相對隆起，而其他區(qū)域均下沉；在隆起幅度方面，靠近起爆點的管廊上覆土層隆起劇烈，而遠離起爆點的管廊上覆土層隆起幅度較小。

圖11 不同時刻地面豎向位移響應分布情況Fig. 11 Distribution of ground vertical displacement at different times

3.2.2 管廊縱向地表豎向位移分布特性

由于淺埋管廊爆炸中地表豎向位移遠大于水平位移，因此著重分析地表豎向位移的響應。圖12為管廊縱向距起爆點水平距離0、2.0、4.0、6.0、8.0 m 處各質(zhì)點豎向位移隨時間的變化曲線。由圖12可知，管廊縱向質(zhì)點豎向位移隨距起爆點水平位移距離的增大而下降明顯。結(jié)合3.1 節(jié)的壓強分析，爆炸過程中管廊縱向地面各質(zhì)點位移響應經(jīng)歷了3 個階段：0～0.05 s，地面質(zhì)點響應壓強快速增加，而豎向位移增加緩慢，其中0 和2.0 m 處質(zhì)點位移向上，4.0、6.0、8.0 m 處質(zhì)點位移向下；0.05～0.39 s，地面質(zhì)點響應壓強度過峰值階段，而質(zhì)點豎向位移增加迅速，靠近起爆點正上方的質(zhì)點由于模型整體沉降作用而呈現(xiàn)先上升后下降的過程；0.39～0.50 s，地面響應壓強逐漸穩(wěn)定至較小值，質(zhì)點位移增加速率減小。

圖13 為管廊縱向以地面上距起爆點最遠處質(zhì)點為參考對象，質(zhì)點相對豎向位移隨起爆點距離的變化曲線。由圖13 可知，管廊縱向相對豎向位移隨距起爆點距離呈較明顯的階段性特征，其中0～2.5 m 處地面向上劇烈隆起，管廊由于甲烷爆炸作用發(fā)生裂損，氣體涌出直接沖擊上覆巖土層；2.5～7.5 m 處為過渡階段，地面隆起程度隨距離的過渡性變化較為明顯；7.5～12.5 m處為平緩階段，地面隆起平緩，這是由于爆炸作用使得管廊整體發(fā)生振動引起地面隆起。整體上，管廊縱向相對豎向位移隨距起爆點水平距離的增大而減小，相對位移峰值出現(xiàn)在起爆點正上方，相對隆起高度為1.3 m。

圖12 管廊縱向質(zhì)點豎向位移隨時間變化曲線Fig. 12 Variation curves of vertical displacement of particles in the longitudinal direction of pipe gallery with time

圖13 管廊縱向質(zhì)點相對豎向位移隨起爆點距離的變化Fig. 13 Relative vertical displacement of particles in the longitudinal direction of pipe gallery with distance of initiation point

3.2.3 管廊橫向地表豎向位移分布特性

圖14 為沿管廊橫向距起爆點水平距離-8.0、-4.0、0、4.0、8.0 m 處質(zhì)點相對豎向位移隨時間的變化曲線。與縱向相似，除了起爆點正上方外，各質(zhì)點豎直位移呈現(xiàn)3 個階段：緩慢增加階段、加速增加階段和減速穩(wěn)定階段。

圖15 為管廊橫向以地面上距起爆點最遠處質(zhì)點為參考對象，質(zhì)點相對豎向位移隨距起爆點距離的變化曲線。管廊橫向相對豎向位移最大值出現(xiàn)在起爆點斜上方，這是由于管廊上覆地面巖土體與管廊黏結(jié)力較小，而巖土體之間的黏結(jié)力使得橫向位移最大值偏移起爆點正上方；各質(zhì)點相對豎向位移隨距起爆點的水平距離增大而減小，其中管廊上覆地面的豎向位移要大于巖層上覆地面，且位移減少速率更??；與各質(zhì)點響應峰值壓強呈現(xiàn)明顯分段性不同，上述散點分布的分段性不明顯，這是由于材料的連續(xù)性；隨著質(zhì)點距起爆點距離的增大，相對豎向位移減少速度逐漸減小。

圖14 管廊橫向質(zhì)點相對豎向位移隨時間變化曲線Fig. 14 Variation curves of relative vertical displacement of particles in the horizontal direction of pipe gallery with time

圖15 管廊橫向質(zhì)點相對豎向位移隨距起爆點距離變化Fig. 15 Rrelative vertical displacement of particles in the horizontal direction of pipe gallery with distance of initiation point

4 結(jié) 論

采用點火增長模型模擬管廊內(nèi)甲烷氣體泄漏爆炸沖擊管廊本體結(jié)構(gòu)和圍巖的過程，研究了地面壓強和位移響應規(guī)律，得出以下主要結(jié)論。

（1）在響應壓強方面：爆炸過程中，管廊及圍巖中出現(xiàn)由接觸面反射和折射產(chǎn)生的次生應力波，這種次生波在管廊橫向方向，振幅隨距起爆點水平距離的增大而增大；而在管廊縱向上產(chǎn)生的次生波振幅較小，且隨距離增大變化較?。还芾葯M縱方向上響應壓強峰值發(fā)生點均偏離起爆點；管廊上覆地面在爆炸作用后出現(xiàn)由地面反射波造成的負向壓強。

（2）在響應位移方面：爆炸作用造成整體地面沉降，但在起爆點中心附近地面隆起，這種隆起由管廊襯砌破裂、氣體直接沖擊巖土體形成的劇烈隆起，以及管廊整體振動形成的輕微隆起兩部分組成，其中圍巖上覆地面位移小于管廊上覆地面；管廊橫向的豎向位移偏移于起爆點正上方。