基于損傷-虛擬張拉裂紋模型的地下爆炸圍巖破壞規律研究*

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

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

球形硐室三維力學性能優越，在軍事工程中應用廣泛。開挖形成的球形腔室，不需要復雜的支護即可作為隱蔽空間進行爆炸實驗。作為爆炸危險物的貯藏點，硐室存在潛在的內爆炸風險，弄清硐室抗爆性能和爆炸后的圍巖損傷規律，對硐室的防護及修復十分重要。

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

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

1 數值模擬方法

1.1 巖石材料的HJC 模型

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

本文研究對象為我國西南地區完整性較好的硬質紅砂巖，通過單軸及常規三軸壓縮實驗、巴西劈裂實驗，得到自然狀態下紅砂巖的基礎力學參數，見表1。

表1 紅砂巖的基本力學參數Table 1 Basic mechanical parameters of red sandstone

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

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

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

圖1 單軸壓縮數值模擬結果Fig. 1 Numerical simulation results of uniaxial compression

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

相應地，為驗證其在動態破壞模擬時的有效性，還進行了SHPB 實驗的數值模擬，建立的實驗裝置幾何模型如圖2 所示。

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

圖2 SHPB 沖擊壓縮實驗的數值模型Fig. 2 Numerical model of the SHPB impact test

圖3 SHPB 沖擊壓縮數值實驗結果Fig. 3 Numerical simulation results of the SHPB impact test

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

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

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

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

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

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

圖5 中，σ為法向應力，為法向剛度，δ為法向相對位移。δ為法向的損傷起始位移，當δ＜δ時，內聚力單元完全彈性，當δ＞δ時，損傷演化開始，δ=δ為極限狀態；σ=δ為法向最大應力；δ為法向的失效位移，當δ＞δ時，虛擬裂紋被激活。

根據圖5(b)，損傷起始后，內聚力單元的本構關系可以表示為：

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

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

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

1.3 模擬方法的可靠性

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

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

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

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

2 地下填實及空腔爆炸的模擬

2.1 算例及建模方法

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

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

數值模擬的建立在顯式動力分析軟件LS-DYNA 中完成，其中巖體部分包括實體單元和節理單元，采用基于Delaunry 三角劃分的Lagrange 網格，空氣域則采用結構化的Euler 網格劃分，空氣域為正方形，并通過關鍵字Initial_Volume_Fraction_Geometry 對炸藥材料空間位置進行定義，以提高網格質量。巖體與空氣重疊部分通過罰函數實現流固耦合，由關鍵字Constrained_Lagrange_In_Solid 控制，Lagrange 網格和Euler 的控制尺寸均為0.2 m，且爆炸產物侵蝕巖體時，不允許網格間出現穿透，巖體外側為無反射邊界。

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

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

選取的TNT 炸藥參數見表3。

表3 TNT 炸藥爆轟產物JWL 參數Table 3 JWL EOS parameters of the TNT detonation product

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

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

在LS-DYNA 中，理想氣體可用附加多項式狀態方程的Null 材料模擬，標準狀態方程為：

式中：=(ρ -ρ/ρ；～為多項式系數，當==0、γ -1、0 時，式(5)與式(4)完全等效。

2.2 爆炸荷載

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

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

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

而對于空氣沖擊波壓力，則可以采用最大擴散速度法進行計算，其表達式為：

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

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

2.3 圍巖分區破壞規律

數值模擬結果見圖9，地下填實爆炸過程中圍巖的損傷在時間和空間上都具有很強的分區特征。

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

根據第1.2 節內容，內聚力單元的失效可以作為裂紋產生和擴展的可靠判據，而實體單元的損傷則可以作為巖石屈服的判據。相應地，對于地下爆炸，兩者可用于裂隙區和破碎區范圍的確定。考慮到裝藥量對爆炸效果的影響，根據霍普金斯比例定律，引入比例半徑來描述地下爆炸的分區破壞規律。比例半徑的定義為：

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

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

圖9 地下填實爆炸的分區破壞規律Fig. 9 Zonal failure law of the underground coupled explosion

2.4 硐室內爆炸的空腔解耦

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

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

3 討 論

對本文損傷-虛擬裂紋模型與傳統FEM 模型的差異與優劣，進行簡單討論。

在圖10(a)所示的耦合裝藥模型的基礎上，通過調整裝藥半徑，可以將裝藥形式轉化為不耦合裝藥，從而減小圍巖損傷，當裝藥半徑調整為0.12 m 時，對應的損傷及虛擬裂紋如圖11(a)～(b)所示，圍巖中破碎區范圍較小，且出現了多條細長裂紋，形成半徑2 m 的裂隙區。

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

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

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

4 結 論

在巖石HJC 材料模型的基礎上，基于結理內聚力單元，提出了一種新型的損傷-虛擬裂紋模型。基于該模型進行了地下填實及空腔爆炸的數值模擬，得到以下結論。

（1）本文的模擬方法不僅能夠很好地模擬巖石的準靜態和動態壓縮過程，而且克服了靜水壓力為拉應力時，HJC 模型破壞形態異常的問題，擴大了HJC 模型的適用范圍，且模型中實體單元的損傷和節理單元的失效可以分別作為巖石壓碎和開裂的判定依據。

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

（3）根據模擬結果，地下填實爆炸具有明顯的分區破壞規律，而隨著硐室尺寸的增大，空腔效應可以減小爆炸荷載對圍巖的損傷作用，對于硬質砂巖硐室，比例半徑為0.52 m/kg時可以實現爆炸荷載的完全解耦，研究結果對地下硐室的防爆抗爆設計有一定的指導作用。