程樹范，葉 陽，曾亞武，高 睿

（武漢大學(xué)土木建筑工程學(xué)院，湖北 武漢 430072）

球形硐室三維力學(xué)性能優(yōu)越，在軍事工程中應(yīng)用廣泛。開挖形成的球形腔室，不需要復(fù)雜的支護(hù)即可作為隱蔽空間進(jìn)行爆炸實(shí)驗(yàn)。作為爆炸危險(xiǎn)物的貯藏點(diǎn)，硐室存在潛在的內(nèi)爆炸風(fēng)險(xiǎn)，弄清硐室抗爆性能和爆炸后的圍巖損傷規(guī)律，對硐室的防護(hù)及修復(fù)十分重要。

與填實(shí)爆炸相比，大型空腔的間隔作用有效規(guī)避了爆炸產(chǎn)物對腔壁的直接沖擊，并削弱了爆炸沖擊波的強(qiáng)度，有利于維持周邊巖體的穩(wěn)定性。李孝蘭對空腔解耦效應(yīng)進(jìn)行了系統(tǒng)總結(jié)，認(rèn)為隨著空腔半徑的增大，爆炸荷載最終能夠?qū)崿F(xiàn)解耦。樓溈濤等在硬巖場地中進(jìn)行了球型空腔的間隔爆炸試驗(yàn)，結(jié)果表明，半徑1.85 m 的球型硐室基本可以實(shí)現(xiàn)240 kg TNT 炸藥的完全解耦；王占江則進(jìn)行了微量裝藥的解耦模擬試驗(yàn)，也驗(yàn)證了空腔解耦的可行性。由于爆炸試驗(yàn)危險(xiǎn)性高，且試驗(yàn)結(jié)果難以被準(zhǔn)確地量測，因此數(shù)值模擬被廣泛用于地下結(jié)構(gòu)的抗爆設(shè)計(jì)，如：Wang 等采用數(shù)值模擬研究了多源爆炸條件下地下硐室的動(dòng)態(tài)響應(yīng)，認(rèn)為應(yīng)加強(qiáng)爆炸源中段的防護(hù)措施；熊益波等基于數(shù)值模擬技術(shù)對抗爆硐室可能存在的高風(fēng)險(xiǎn)區(qū)進(jìn)行了預(yù)測，并制定了有針對性的加強(qiáng)方案。對于爆炸問題，在連續(xù)介質(zhì)理論的基礎(chǔ)上，針對巖體大變形和高應(yīng)變率工況，研究人員有針對性地提出并改進(jìn)了材料模型，使得模擬結(jié)果更加可靠，然而在模擬脆性開裂時(shí)仍存在明顯的不足。而離散元方法雖然可以較好地模擬裂紋的形成與擴(kuò)展，但并不適合大尺度的模擬。目前，部分研究基于“生死”單元，在有限元模型中通過單元失效模擬裂紋，取得了較好的模擬效果，但單元?jiǎng)h除導(dǎo)致的質(zhì)量不守恒和能量異常損失仍會(huì)降低數(shù)值模擬的可靠性，特別是當(dāng)網(wǎng)格尺寸較大或破碎程度較嚴(yán)重時(shí)，該方法的模擬結(jié)果將嚴(yán)重失真。

為克服連續(xù)介質(zhì)方法的不足，本文中在巖石HJC 模型的基礎(chǔ)上，提出一種新型的損傷-虛擬裂紋模型，首先討論該模型的適用性和尺寸效應(yīng)，并基于該模型對硐室內(nèi)爆炸引起的圍巖損傷規(guī)律進(jìn)行分析，最后通過與有限元方法的對比，討論該模型及模擬方法的優(yōu)越性。

1 數(shù)值模擬方法

1.1 巖石材料的HJC 模型

數(shù)值模擬過程中，材料模型的選擇是否合理直接決定了數(shù)值模擬的可靠性，Holmquist 等在金屬材料Johnson-Cook 模型基礎(chǔ)上提出的HJC 模型綜合考慮了靜水壓力、應(yīng)變率和損傷對材料強(qiáng)度及力學(xué)性能的影響，是模擬脆性材料動(dòng)態(tài)破壞的一種較理想的模型。HJC 模型有19 個(gè)獨(dú)立的模型參數(shù)，分別為基本力學(xué)參數(shù)、、、ρ，極限面參數(shù)、、、，壓力參數(shù)、、、、、、，應(yīng)變率參數(shù)和損傷參數(shù)ε、、。HJC 模型雖然包含較多參數(shù)，但其物理意義明確，標(biāo)定方法也已經(jīng)較為成熟。

本文研究對象為我國西南地區(qū)完整性較好的硬質(zhì)紅砂巖，通過單軸及常規(guī)三軸壓縮實(shí)驗(yàn)、巴西劈裂實(shí)驗(yàn)，得到自然狀態(tài)下紅砂巖的基礎(chǔ)力學(xué)參數(shù)，見表1。

表1 紅砂巖的基本力學(xué)參數(shù)Table 1 Basic mechanical parameters of red sandstone

表1 中單軸壓縮實(shí)驗(yàn)對應(yīng)的加載應(yīng)變率為10s，為標(biāo)定應(yīng)變率參數(shù)，另外還進(jìn)行了3 組不同應(yīng)變率的單軸壓縮實(shí)驗(yàn)和1 組分離式霍普金森壓桿(split Hopkinson pressure bar，SHPB)實(shí)驗(yàn)，得到應(yīng)變率為5×10、2×10、10和87.3 s時(shí)硬質(zhì)砂巖的動(dòng)態(tài)強(qiáng)度分別為79.18、83.57、85.13 和100.65 MPa。最后根據(jù)文獻(xiàn)[18-19]中所建議的方法，得到HJC 模型參數(shù)，見表2。

表2 紅砂巖HJC 模型參數(shù)Table 2 Parameters of the HJC model of red sandstone

采用表2 中的材料參數(shù)，基于有限元方法（FEM）在LS-DYNA 軟件平臺(tái)上進(jìn)行了單軸壓縮的數(shù)值模擬，以驗(yàn)證所標(biāo)定參數(shù)的可靠性?？紤]到數(shù)值計(jì)算的效率，根據(jù)對稱性，將圓柱體試件簡化為圖1(a)所示的二維模型進(jìn)行分析，模型尺寸為50 mm×100 mm，采用非結(jié)構(gòu)化的Delaunry 網(wǎng)格劃分，網(wǎng)格控制尺寸為2 mm，壓力板與試件間采用基于罰函數(shù)的雙向接觸，接觸摩擦因數(shù)為0.1。模擬得到的單軸壓縮應(yīng)力-應(yīng)變曲線和破壞后的形態(tài)如圖1(b)所示。

圖1 單軸壓縮數(shù)值模擬結(jié)果Fig. 1 Numerical simulation results of uniaxial compression

由圖1(b)可知，單軸壓縮數(shù)值模擬得到的應(yīng)力-應(yīng)變曲線與實(shí)驗(yàn)值一致性較好。雖然FEM 模型的最終破壞形態(tài)與實(shí)際狀態(tài)存在一些差異，但均以剪切破壞為主，且模擬得到的破壞傾角與主裂紋傾角十分接近，對應(yīng)的破壞機(jī)理是一致的。可見經(jīng)過參數(shù)標(biāo)定的HJC 模型在模擬巖石準(zhǔn)靜態(tài)單軸壓縮時(shí)的結(jié)果是可靠的。

相應(yīng)地，為驗(yàn)證其在動(dòng)態(tài)破壞模擬時(shí)的有效性，還進(jìn)行了SHPB 實(shí)驗(yàn)的數(shù)值模擬，建立的實(shí)驗(yàn)裝置幾何模型如圖2 所示。

圖2 中的入射桿和透射桿均為低合金鋼材質(zhì)，彈性模量為210 GPa，泊松比為0.27。入射桿、透射桿的長度分別為2.5、1.5 m，直徑為50 mm；巖石試件高度和直徑均為50 mm，桿件和試件的網(wǎng)格尺寸分別為2、5 mm。對于入射應(yīng)力波的加載，采取時(shí)程荷載的方式來實(shí)現(xiàn)，即將實(shí)驗(yàn)測得的入射波直接施加于入射桿，因此數(shù)值模型中并不包括子彈。實(shí)驗(yàn)與數(shù)值模擬得到的透射波、入（反）射波以及由三波法得到的應(yīng)力-應(yīng)變曲線如圖3 所示。

圖2 SHPB 沖擊壓縮實(shí)驗(yàn)的數(shù)值模型Fig. 2 Numerical model of the SHPB impact test

圖3 SHPB 沖擊壓縮數(shù)值實(shí)驗(yàn)結(jié)果Fig. 3 Numerical simulation results of the SHPB impact test

由于采用了時(shí)程荷載的方式施加入射波，圖3(a)中入射波段完全重合，而模擬得到的反射波和透射波也具有較高的吻合度，數(shù)值模擬得到的加載平均應(yīng)變率為90.36 s，動(dòng)態(tài)抗壓強(qiáng)度為103.22 MPa，相應(yīng)的實(shí)驗(yàn)值分別為87.3 s和100.65 MPa，誤差均小于5%，說明經(jīng)標(biāo)定的HJC 模型在模擬動(dòng)態(tài)破壞時(shí)也具有較為理想的效果。

然而在本構(gòu)關(guān)系方面，相較于壓縮破壞時(shí)的分段考慮，HJC 模型在模擬巖石受拉破壞時(shí)采用的理想彈塑性本構(gòu)并不適合描述巖石的脆性斷裂。采用表2 中的參數(shù)，基于FEM 模擬了巴西劈裂實(shí)驗(yàn)，模擬方法與單軸壓縮試驗(yàn)相同。得到的破壞形態(tài)及應(yīng)力-應(yīng)變曲線如圖4 所示，與室內(nèi)實(shí)驗(yàn)有較大的差異，試件最終表現(xiàn)為與單軸壓縮相似的斜裂紋破壞。這是由于HJC 模型僅在受壓時(shí)考慮了損傷對材料力學(xué)性能的影響，而損傷變量的定義又包含了張拉產(chǎn)生的塑性應(yīng)變，在靜水壓力為拉的情況下，張拉產(chǎn)生的損傷會(huì)導(dǎo)致材料抗壓強(qiáng)度快速下降，產(chǎn)生圖4 所示的異常破壞形式，這是HJC 模型的一個(gè)重大不足。

圖4 采用有限元方法進(jìn)行的巴西劈裂數(shù)值模擬Fig. 4 Numerical test of Brazilian splitting test based on FEM

1.2 考慮張拉失效的虛擬裂紋模型

為克服HJC 模型在模擬張拉破壞時(shí)的不足，考慮引入虛擬裂紋來模擬張拉破壞，而不考慮材料本身的張拉屈服，即采用如圖5(a)所示的節(jié)理單元連接代替FEM 中所采用的實(shí)體單元共節(jié)點(diǎn)連接。建模過程中，首先將單元離散化，此時(shí)節(jié)點(diǎn)在空間上重合，但各自獨(dú)立，其后在相鄰實(shí)體單元間插入一無厚度節(jié)理單元，通過節(jié)理單元與實(shí)體單元間的共節(jié)點(diǎn)連接，最后將離散單元重新組合為整體。節(jié)理單元采用張拉失效的雙線性內(nèi)聚力模型，內(nèi)聚力單元受力后可以在法向產(chǎn)生獨(dú)立的位移，而在切向不會(huì)發(fā)生變形，其法向應(yīng)力-位移曲線如圖5(b)所示，包括無損彈性和損傷軟化兩個(gè)階段。

圖5 張拉失效的雙線性內(nèi)聚力模型Fig. 5 Cohesive model with a tension failure rule

圖5 中，σ為法向應(yīng)力，為法向剛度，δ為法向相對位移。δ為法向的損傷起始位移，當(dāng)δ＜δ時(shí)，內(nèi)聚力單元完全彈性，當(dāng)δ＞δ時(shí)，損傷演化開始，δ=δ為極限狀態(tài)；σ=δ為法向最大應(yīng)力；δ為法向的失效位移，當(dāng)δ＞δ時(shí)，虛擬裂紋被激活。

根據(jù)圖5(b)，損傷起始后，內(nèi)聚力單元的本構(gòu)關(guān)系可以表示為：

式中：為法向的損傷變量，可定義為：

由于節(jié)理單元在切向不會(huì)出現(xiàn)位移，剪切變形需要通過實(shí)體單元來完成，因此，在壓縮狀態(tài)下，本文模型將退化為FEM 模型。插入節(jié)理單元后，將引入3 個(gè)細(xì)觀材料參數(shù)、δ和δ，其取值無法通過實(shí)驗(yàn)直接得到，需要通過巴西劈裂的數(shù)值模擬反演得到，即要求模擬得到的荷載-位移曲線以及最終的破壞模式均與實(shí)驗(yàn)一致。最終標(biāo)定的細(xì)觀參數(shù)為：=4.2 GPa/mm，δ=0.012 mm，δ=0.025 mm，對應(yīng)的荷載-位移曲線如圖6(a)所示，無論是應(yīng)力-應(yīng)變曲線，還是最終的破壞狀態(tài)，均與實(shí)驗(yàn)高度吻合。根據(jù)數(shù)值模擬結(jié)果，內(nèi)聚力單元對單軸壓縮過程的模擬結(jié)果影響較小，插入內(nèi)聚力單元后的應(yīng)力-應(yīng)變曲線和破壞形態(tài)如圖6(b)所示。

圖6 巴西劈裂及單軸壓縮數(shù)值模擬Fig. 6 Numerical simulation of the Brazilian splitting and uniaxial compression based on the proposed method

1.3 模擬方法的可靠性

虛擬裂紋的引入，雖然解決了HJC 模型在模擬低靜水壓力時(shí)的不足，但虛擬裂紋的插入不可避免地會(huì)導(dǎo)致數(shù)值模型出現(xiàn)網(wǎng)格依賴性和尺寸效應(yīng)。且插入節(jié)理單元后，數(shù)值模型的計(jì)算量將顯著增加，對于大型工程尺度的模擬計(jì)算，難以實(shí)現(xiàn)毫米級的網(wǎng)格劃分，增大網(wǎng)格尺寸后，標(biāo)定的模型參數(shù)是否適用，還需要進(jìn)行專門的討論。

HJC 模型的材料參數(shù)本身與單元尺寸無關(guān)，在模型被放大后，理論上不需要重新標(biāo)定，但內(nèi)聚力單元的參數(shù)與網(wǎng)格尺寸密切相關(guān)，對于張拉破壞，由于破壞面較為穩(wěn)定，可以引入一個(gè)比例因子來考慮尺寸效應(yīng)，即在保持破壞應(yīng)變不變的情況下，調(diào)整模量和位移來控制尺寸效應(yīng)。當(dāng)網(wǎng)格尺寸擴(kuò)大倍后，對應(yīng)的法向剛度參數(shù)保持不變，法向位移參數(shù)δ和δ則同步擴(kuò)大倍，這樣處理即可保持變形模量和破壞應(yīng)變不變。為驗(yàn)證這一方法的可靠性，在保持網(wǎng)格形式一致的前提下，分別建立50 mm×100 mm、50 cm×100 cm、50 dm×100 dm 以及50 m×100 m 等4 個(gè)尺度的單軸壓縮模型和單軸拉伸模型，對應(yīng)的比例因子分別為1、10、100 和1 000。將單軸加載速率依次設(shè)置為2.5 mm/s、2.5 cm/s、2.5 dm/s 和2.5 m/s，以保證加載應(yīng)變率一致，數(shù)值模擬得到的應(yīng)力-應(yīng)變曲線如圖7 所示。

圖7 不同尺度的單軸壓縮和單軸拉伸應(yīng)力-應(yīng)變曲線Fig. 7 Stress-strain curves of uniaxial compression and uniaxial tension at different scales

由于本文模型的剪切破壞受節(jié)理單元的影響較小，在保持網(wǎng)格劃分形式不變的基礎(chǔ)上，隨著網(wǎng)格尺寸的增大，材料強(qiáng)度和變形模量僅出現(xiàn)了小幅度的下降。在圖7(a)中，當(dāng)網(wǎng)格尺寸為2 cm 和2 dm 時(shí)，模型的強(qiáng)度為2 mm 網(wǎng)格時(shí)的94.55%和91.75%，減小幅度未超過10%，因此采用厘米和分米級的網(wǎng)格進(jìn)行工程模擬也具有較高的可行性，且10%以內(nèi)的誤差可以作為安全儲(chǔ)備，模擬結(jié)果偏于安全。當(dāng)網(wǎng)格尺寸達(dá)到米的量級時(shí)，誤差仍能保持在20%以內(nèi)，模擬結(jié)果尚可用于定性分析和粗略估算。根據(jù)圖7(b)，當(dāng)網(wǎng)格尺度由毫米量級增加至米量級后，相應(yīng)的抗拉強(qiáng)度增加僅為4.2%，說明本文模型通過比例因子法較好地解決了虛擬裂紋帶來的網(wǎng)格尺寸效應(yīng)，可適用于多種尺度網(wǎng)格。由于本文方法同時(shí)考慮了巖石材料在高靜水壓力狀態(tài)下的應(yīng)變率效應(yīng)、塑性損傷和低靜水壓力狀態(tài)下的脆性斷裂，因此適合于大尺度的地下爆炸的模擬。

2 地下填實(shí)及空腔爆炸的模擬

2.1 算例及建模方法

根據(jù)文獻(xiàn)[6]，硬巖硐室的完全解耦（圍巖對爆炸荷載的響應(yīng)為完全彈性）時(shí)的比例半徑約為0.35 m/kg，那么當(dāng)TNT 炸藥裝藥半徑=0.5 m（裝藥量863 kg）時(shí)，解耦半徑約為3.3 m。據(jù)此分別建立填實(shí)爆炸模型，和空腔半徑=1、2、3.5 和5 m 的4 個(gè)空腔爆炸模型。根據(jù)試算，巖體部分模型尺寸取為20 m×20 m，空氣域范圍為6 m×6 m，即可實(shí)現(xiàn)爆炸荷載的有效耦合且虛擬裂紋及損傷均不會(huì)超出模型區(qū)域，以=2 m為例，建立的數(shù)值模型如圖8 所示。

圖8 空腔爆炸數(shù)值模型(R=2 m)Fig. 8 Numerical model of cavity explosion in LS-DYNA (R=2 m)

數(shù)值模擬的建立在顯式動(dòng)力分析軟件LS-DYNA 中完成，其中巖體部分包括實(shí)體單元和節(jié)理單元，采用基于Delaunry 三角劃分的Lagrange 網(wǎng)格，空氣域則采用結(jié)構(gòu)化的Euler 網(wǎng)格劃分，空氣域?yàn)檎叫?，并通過關(guān)鍵字Initial_Volume_Fraction_Geometry 對炸藥材料空間位置進(jìn)行定義，以提高網(wǎng)格質(zhì)量。巖體與空氣重疊部分通過罰函數(shù)實(shí)現(xiàn)流固耦合，由關(guān)鍵字Constrained_Lagrange_In_Solid 控制，Lagrange 網(wǎng)格和Euler 的控制尺寸均為0.2 m，且爆炸產(chǎn)物侵蝕巖體時(shí)，不允許網(wǎng)格間出現(xiàn)穿透，巖體外側(cè)為無反射邊界。

LS-DYNA 軟件自身擁有豐富的材料庫，其中High_Explosive_Burn 材料被廣泛用于凝聚相炸藥的模擬，以TNT 炸藥為例，其爆速=6 930 m/s，爆壓=21 GPa。爆炸產(chǎn)物的擴(kuò)容過程則通過附加JWL（Jones-Wilkins-Lee）狀態(tài)方程來實(shí)現(xiàn)，其表達(dá)式為：

式中：、為線性爆炸系數(shù)，具有壓強(qiáng)的量綱；ω、和為無量綱爆炸系數(shù)；為相對體積，即壓強(qiáng)為時(shí)氣體體積與初始體積之比；為炸藥單位體積爆轟能量。

選取的TNT 炸藥參數(shù)見表3。

表3 TNT 炸藥爆轟產(chǎn)物JWL 參數(shù)Table 3 JWL EOS parameters of the TNT detonation product

炸藥為中心起爆，爆炸過程空氣可視為理想氣體，其狀態(tài)方程為：

式中：為空氣壓力；ρ/ρ為相對密度，其中ρ 為壓縮后空氣的密度，ρ為空氣初始密度，取為1.29 kg/m；γ 為理想氣體比熱容，取為1.4；為單位體積內(nèi)能，取為250 kJ/m。

在LS-DYNA 中，理想氣體可用附加多項(xiàng)式狀態(tài)方程的Null 材料模擬，標(biāo)準(zhǔn)狀態(tài)方程為：

式中：=(ρ -ρ/ρ；～為多項(xiàng)式系數(shù)，當(dāng)==0、γ -1、0 時(shí)，式(5)與式(4)完全等效。

2.2 爆炸荷載

對于地下填實(shí)爆炸，爆炸荷載以爆炸產(chǎn)物壓力的形式直接作用于巖體，其峰值壓力可達(dá)10 GPa 數(shù)量級；而對于空腔爆炸，由于空氣的間隔作用，爆炸荷載多以空氣沖擊波的形式作用于巖體，數(shù)量級相對較低。填實(shí)和空腔爆炸在荷載形式上有較大區(qū)別，荷載的計(jì)算方法也有所差異。

對于填實(shí)爆炸，爆炸產(chǎn)物壓力可以采用等熵膨脹法進(jìn)行計(jì)算，亨利奇給出的耦合裝藥最大壓力的計(jì)算方法為：

式中：為作用于巖體的最大壓力；ρ為炸藥密度，取為1 650 kg/m（TNT）；為炸藥爆轟速度，取為6 930 m/s；γ 為等熵指數(shù)，取為3。

而對于空氣沖擊波壓力，則可以采用最大擴(kuò)散速度法進(jìn)行計(jì)算，其表達(dá)式為：

式中：ρ為空氣密度；為巖石內(nèi)沖擊波速，當(dāng)沖擊波強(qiáng)度較小時(shí)，可取為剪切波速1 264 m/s；為空氣絕熱指數(shù)，對于雙原子分子可取為1.41；為沖擊波壓力增大系數(shù)，取值為0～20。

以=0.5， 5 m 的模型為例，由式(6) 計(jì)算得到的填實(shí)爆炸對應(yīng)的孔壁最大壓力為9.91 GPa；由式(7)計(jì)算得到的空腔爆炸最大壓力范圍為0～34.2 MPa。而本文數(shù)值模擬監(jiān)測到的孔壁最大壓力分別為8.45 GPa 和10.3 MPa，均在理論值附近，說明ALE 算法能夠較為準(zhǔn)確地模擬爆炸荷載作用。

2.3 圍巖分區(qū)破壞規(guī)律

數(shù)值模擬結(jié)果見圖9，地下填實(shí)爆炸過程中圍巖的損傷在時(shí)間和空間上都具有很強(qiáng)的分區(qū)特征。

在時(shí)間上可以分為3 個(gè)階段：第1 階段為孔壁的外擴(kuò)，此階段內(nèi)爆炸氣體直接沖擊巖體，產(chǎn)生極大的壓力，遠(yuǎn)超巖石的屈服強(qiáng)度，巖石表現(xiàn)出一種類似于流體的力學(xué)性質(zhì)，巖體向外移動(dòng)形成爆炸空腔，并激發(fā)壓縮波；第2 階段為破碎區(qū)的形成，此階段內(nèi)壓縮波強(qiáng)度大于巖石的屈服強(qiáng)度，巖體出現(xiàn)大的塑性損傷，巖體完全破碎；第3 階段為裂隙區(qū)的形成和擴(kuò)展，此階段內(nèi)壓縮波壓力小于巖體的抗壓強(qiáng)度，但受到泊松效應(yīng)影響，環(huán)向的拉應(yīng)力仍可誘發(fā)張拉裂紋的擴(kuò)展。隨著壓縮波壓力下降至不足以誘發(fā)裂紋繼續(xù)擴(kuò)展后，爆炸過程結(jié)束。相應(yīng)地，根據(jù)破壞機(jī)制，在空間上則會(huì)形成擴(kuò)容空腔（見圖9(a)）、壓縮破壞區(qū)和張拉破壞區(qū)（見圖9(c)）。

根據(jù)第1.2 節(jié)內(nèi)容，內(nèi)聚力單元的失效可以作為裂紋產(chǎn)生和擴(kuò)展的可靠判據(jù)，而實(shí)體單元的損傷則可以作為巖石屈服的判據(jù)。相應(yīng)地，對于地下爆炸，兩者可用于裂隙區(qū)和破碎區(qū)范圍的確定?？紤]到裝藥量對爆炸效果的影響，根據(jù)霍普金斯比例定律，引入比例半徑來描述地下爆炸的分區(qū)破壞規(guī)律。比例半徑的定義為：

式中：為爆心距；為裝藥質(zhì)量，對于本文模型，可取=863 kg。

如圖9(c)所示，耦合裝藥時(shí)破碎區(qū)平均半徑為2.0 m，是炮孔半徑的4 倍；考慮裝藥尺寸后的比例半徑為0.26 m/kg；而裂隙區(qū)半徑為5.0 m，是炮孔半徑的10 倍，對應(yīng)的比例半徑為0.65 m/kg。

圖9 地下填實(shí)爆炸的分區(qū)破壞規(guī)律Fig. 9 Zonal failure law of the underground coupled explosion

2.4 硐室內(nèi)爆炸的空腔解耦

數(shù)值模擬得到的空腔爆炸模型，其最終的破壞形式如圖10 所示。

對于空腔爆炸，破碎區(qū)的厚度隨著空腔尺寸的增大而減小，說明空腔效應(yīng)可以減小作用于孔壁上的爆炸壓力，從而提高硐室的抗爆性能。當(dāng)硐室半徑達(dá)到預(yù)設(shè)的解耦尺寸（=3.5 m）時(shí)，圍巖表面僅出現(xiàn)零星的實(shí)體單元損傷，而壁面附近卻分布有密集的短裂紋。分析認(rèn)為，當(dāng)爆炸沖擊波在硐室表面反射時(shí)會(huì)產(chǎn)生反射，反射波與后續(xù)的沖擊波疊加使爆炸壓力被加強(qiáng)，并轉(zhuǎn)化為張拉應(yīng)力波，由于巖體的抗拉強(qiáng)度遠(yuǎn)小于抗壓強(qiáng)度，張拉作用仍能引起表層的局部破壞，這一現(xiàn)象在文獻(xiàn)[10]中也被觀測和驗(yàn)證，這部分由反射波疊加形成的裂紋雖然十分密集但沒有擴(kuò)展的趨勢，因此仍將其歸于裂隙區(qū)。當(dāng)硐室半徑達(dá)到5 m 后，裂隙區(qū)完全消失，硐室處于完全彈性狀態(tài)，據(jù)此認(rèn)為：對于硬質(zhì)砂巖硐室，比例半徑為0.52 m/kg時(shí)，可以實(shí)現(xiàn)爆炸荷載的完全解耦。

3 討 論

對本文損傷-虛擬裂紋模型與傳統(tǒng)FEM 模型的差異與優(yōu)劣，進(jìn)行簡單討論。

在圖10(a)所示的耦合裝藥模型的基礎(chǔ)上，通過調(diào)整裝藥半徑，可以將裝藥形式轉(zhuǎn)化為不耦合裝藥，從而減小圍巖損傷，當(dāng)裝藥半徑調(diào)整為0.12 m 時(shí)，對應(yīng)的損傷及虛擬裂紋如圖11(a)～(b)所示，圍巖中破碎區(qū)范圍較小，且出現(xiàn)了多條細(xì)長裂紋，形成半徑2 m 的裂隙區(qū)。

圖10 空腔爆炸模型最終的破壞形態(tài)Fig. 10 Final failure characteristics of the cavity explosion models

圖11(c)為相同條件（除無黏聚力單元外，其他條件均與本文模型相同）的有限元方法模擬的損傷云圖。對比圖11(b)～(c)可知：FEM 模型計(jì)算得到的損傷具有很強(qiáng)的對稱性，損傷云圖近似呈環(huán)狀，而本文的損傷-虛擬裂紋模型計(jì)算得到的損傷演化受裂紋擴(kuò)展的影響較明顯，裂紋延伸處的損傷更加嚴(yán)重，對應(yīng)云圖呈現(xiàn)放射狀，體現(xiàn)了裂紋對損傷演化過程的導(dǎo)向作用，虛擬裂紋的隨機(jī)性即可反映出巖石的各向異性特征。由于節(jié)理單元失效消耗了一部分能量，有限元總損傷面積和損傷程度相對較小，但由于裂紋分布不均勻，實(shí)際的損傷演化范圍較有限元模型略大。綜上所述，損傷-虛擬裂紋模型能夠考慮到巖體裂紋對損傷分布的導(dǎo)向作用，綜合考慮損傷及裂紋效應(yīng)后的模擬結(jié)果在理論上應(yīng)該更加真實(shí)。

圖11 本文模擬方法與有限元方法的對比Fig. 11 Comparison between the method in this paper and the FEM

4 結(jié) 論

在巖石HJC 材料模型的基礎(chǔ)上，基于結(jié)理內(nèi)聚力單元，提出了一種新型的損傷-虛擬裂紋模型?；谠撃Ｐ瓦M(jìn)行了地下填實(shí)及空腔爆炸的數(shù)值模擬，得到以下結(jié)論。

（1）本文的模擬方法不僅能夠很好地模擬巖石的準(zhǔn)靜態(tài)和動(dòng)態(tài)壓縮過程，而且克服了靜水壓力為拉應(yīng)力時(shí)，HJC 模型破壞形態(tài)異常的問題，擴(kuò)大了HJC 模型的適用范圍，且模型中實(shí)體單元的損傷和節(jié)理單元的失效可以分別作為巖石壓碎和開裂的判定依據(jù)。

（2）通過引入比例因子后，本文模型可以避免不同網(wǎng)格尺寸下對模型參數(shù)的重復(fù)標(biāo)定，使得本文模型可用于工程尺度的分析。同時(shí)本文模型能夠考慮到巖體裂紋對損傷分布的導(dǎo)向作用，綜合考慮損傷及裂紋效應(yīng)后的模擬結(jié)果在理論上也將更加真實(shí)。

（3）根據(jù)模擬結(jié)果，地下填實(shí)爆炸具有明顯的分區(qū)破壞規(guī)律，而隨著硐室尺寸的增大，空腔效應(yīng)可以減小爆炸荷載對圍巖的損傷作用，對于硬質(zhì)砂巖硐室，比例半徑為0.52 m/kg時(shí)可以實(shí)現(xiàn)爆炸荷載的完全解耦，研究結(jié)果對地下硐室的防爆抗爆設(shè)計(jì)有一定的指導(dǎo)作用。