李曉靜，韓雪瑞，張向陽，劉宇軒

(1.山東建筑大學(xué)土木工程學(xué)院，山東 濟(jì)南 250101；2.山東軌道交通勘察設(shè)計(jì)院有限公司，山東 濟(jì)南 250101)

0 引言

鉆爆法因其性價(jià)比高、適應(yīng)范圍廣，已廣泛地應(yīng)用于采礦、隧道建設(shè)中[1－2]。 爆炸產(chǎn)生的能量釋放不可避免地會(huì)導(dǎo)致圍巖破裂，盡管對(duì)爆炸引起的巖石開裂已開展了深入的研究，但是由于復(fù)雜的化學(xué)反應(yīng)和昂貴的現(xiàn)場實(shí)驗(yàn)材料，其爆炸作用很難通過測量手段獲取，而且隨著開采深度的增加，地應(yīng)力對(duì)爆炸引起的斷裂影響不容忽視。

隨著計(jì)算機(jī)運(yùn)算能力的提升，數(shù)值軟件為研究巖石的爆破提供了一種新的途徑。 有限元法(Finite Element Method， FEM)[3]、 有限差分法(Finite Difference Method，F(xiàn)DM)[4]和離散元法(Discrete Element Method，DEM)[5]等數(shù)值模擬方法，均已廣泛應(yīng)用于巖體中爆破的研究。 Banadaki 等[6]通過實(shí)驗(yàn)室爆破試驗(yàn)研究了沖擊波產(chǎn)生的鉆孔周邊的3 個(gè)不同裂紋區(qū)域，測量各區(qū)域力學(xué)性能和爆破壓力峰值及所引起的斷裂模式，并采用有限元程序AUTODYN 校準(zhǔn)了壓力和動(dòng)態(tài)斷裂，由數(shù)值分析所得到的斷裂模式與試驗(yàn)結(jié)果一致。 甯尤軍等[7]開發(fā)了一種不連續(xù)變形分析的方法(Discontinuous Deformation Analysis ，DDA)模擬巖石破壞問題，并有效地預(yù)測了裂縫的萌生、擴(kuò)展和合并過程。 謝理想等[8]通過有限元程序LS－DYNA 數(shù)值模擬了地應(yīng)力作用下的巖體爆破行，表明壓力荷載作用會(huì)影響巖體損傷范圍，隨著地應(yīng)力增大，掏槽爆破巖體的損傷范圍會(huì)變小。 袁增森等[9]通過顆粒流程序(Particle Flow Code，PFC)研究了不耦合裝藥對(duì)花崗巖爆破損傷程度的影響，得出了隨著裝藥不耦合系數(shù)的增大，花崗巖爆破損傷程度先增強(qiáng)后減弱的結(jié)果。 謝冰等[10]運(yùn)用有限元軟件AUTODYN 和離散元軟件(Universal Distinct Element Code，UDEC)相結(jié)合的方法研究了節(jié)理幾何特征對(duì)爆破的影響，并指出節(jié)理間距對(duì)預(yù)裂爆破有較大影響。 在這些方法中，雖然顯式動(dòng)態(tài)有限元軟件可用于巖石爆破建模，但其連續(xù)性假設(shè)會(huì)使有限元難以模擬巖石中裂縫的萌生和巖石破碎，特別是對(duì)于爆炸載荷下大規(guī)模復(fù)雜裂縫的模擬。 與有限元法相比，離散元法在描述巖石微觀結(jié)構(gòu)和模擬巖石破裂過程方面具有明顯的優(yōu)勢。

文章以花崗巖為研究對(duì)象，采用LS－DYNA 和PFC2D 相結(jié)合的數(shù)值分析方法，研究了爆炸沖擊波對(duì)不同應(yīng)力狀態(tài)下花崗巖中裂紋產(chǎn)生和擴(kuò)展的作用機(jī)理，為深部巖體爆破工程提供理論支持。

1 LS－DYNA 和PFC2D 結(jié)合數(shù)值分析方法

1.1 模擬原理

離散元相對(duì)于有限元對(duì)于解決巖石材料的破碎和裂紋的擴(kuò)展具有獨(dú)特的優(yōu)勢，因此，PFC2D 已成為分析巖體裂隙特性的強(qiáng)大工具。 但是，由于離散元程序PFC2D 不能模擬爆炸過程，通常爆炸波的輸入要由其他方法得到，以前的學(xué)者通常采用三角形脈沖或者高斯函數(shù)曲線模擬爆炸應(yīng)力波[11－12]，與真實(shí)的爆炸沖擊波相比有較大的誤差。 文章通過有限元軟件LS－DYNA 模擬爆炸過程，通過監(jiān)測孔壁附近單元，記錄爆炸沖擊速度，將速度轉(zhuǎn)化為PFC 可兼容的數(shù)據(jù)格式，同時(shí)將速度賦予孔壁顆粒，以模擬爆炸應(yīng)力波作用，觀察巖體中的破碎情況及裂縫演化。 簡而言之，沖擊波是通過賦予粒子徑向速度間接施加的。

在PFC2D 程序中，顆粒間的相互作用遵循力－位移定律和牛頓第二定律，巖石材料在PFC2D 中以剛性圓盤組合表示，圓盤之間通過接觸模型粘結(jié)在一起，PB 接觸模型[13]在顆粒之間可以傳遞力和力矩，所以人們普遍認(rèn)為PB 接觸模型是描述巖石微觀特性的最優(yōu)模型。 因此，文章所有試件的顆粒接觸模型均選擇PB 接觸模型。 PB 接觸模型基本部件和力學(xué)性能如圖1 所示。

圖1 PB 接觸模型的基本部件和原理示意圖

PB 接觸模型由顆粒、線性接觸和平行鍵3 個(gè)基本組件組成，可以承受法向力Fn、剪力Fs和彎矩M，其力學(xué)行為可認(rèn)為類似于梁的。 因此，最大拉伸應(yīng)力σt，max和最大剪切應(yīng)力τmax分別由式(1)和(2)表示為

式中A 為粘結(jié)鍵橫截面的面積，mm2；I 為粘結(jié)鍵的慣性矩，mm4；為接觸顆粒中半徑較小值，mm。

如果σt，max超過法向強(qiáng)度或τmax超過剪切強(qiáng)度時(shí)，平行鍵破壞，顆粒間的力和力矩消失，其中微裂紋代表平行鍵斷裂的斷裂，斷裂的平行鍵處代表顆粒之間不再傳遞力和力矩，微裂紋的數(shù)量及位置可以反映巖體的破壞程度。

1.2 計(jì)算模型的建立

以文獻(xiàn)[6]中的室內(nèi)試驗(yàn)作為對(duì)比，制備直徑為144 mm、高度為150 mm 的Barre 花崗巖圓柱形樣品，用于實(shí)驗(yàn)室爆破實(shí)驗(yàn)，如圖2 所示。 在試樣中間鉆取一個(gè)直徑為6.45 mm 的孔作為爆破孔，將厚度為0.6 mm 的銅管緊密安裝在爆孔中，采用直徑為1.65 mm 的太恩(PETN)炸藥及DYNO 公司的導(dǎo)爆索，將空氣作為耦合介質(zhì)進(jìn)行測試，銅管的作用是防止爆生氣體進(jìn)入裂縫，影響對(duì)沖擊波爆破作用的分析，而銅在施加的沖擊載荷下容易變形，并隨鉆孔膨脹而不撕裂，從而防止氣體滲透到產(chǎn)生的裂縫中。

圖2 文獻(xiàn)[6]中爆破試驗(yàn)樣本圖

在軟件分析中精確還原試驗(yàn)的材料性質(zhì)，建模過程的重點(diǎn)在于離散元模型的微觀參數(shù)標(biāo)定和爆破過程的記錄方法。 使用離散元軟件進(jìn)行數(shù)值模擬的一個(gè)關(guān)鍵問題是建立顆粒的微觀參數(shù)與巖石宏觀力學(xué)參數(shù)之間的關(guān)系。 微觀參數(shù)包括法向和剪切剛度、拉伸強(qiáng)度、內(nèi)聚力等，使模型能夠反映真實(shí)巖石試樣的力學(xué)行為。 利用單軸壓縮試驗(yàn)和巴西劈裂試驗(yàn)校準(zhǔn)文獻(xiàn)[6]中的宏觀巖石參數(shù)，通過對(duì)微觀參數(shù)迭代調(diào)整[14]，使模型宏觀參數(shù)與實(shí)際試樣參數(shù)相匹配。 實(shí)際試樣巴雷花崗巖宏觀力學(xué)性能參數(shù)包括密度為2.66 kg/m3、泊松比為0.16，而體積模量、剪切模量分別為25.7 和21.9 GPa，單軸抗拉強(qiáng)度、抗壓強(qiáng)度分別為7.3、161.5 MPa。

巖石的單軸壓縮和巴西劈裂強(qiáng)度模擬試樣如圖3 所示，采用由8 165 個(gè)顆粒組成的寬為150 mm、高為300 mm 的矩形樣品用于模擬單軸壓縮試驗(yàn)，使用895 個(gè)顆粒組成的半徑為75 mm 的圓形樣品用于模擬巴西劈裂試驗(yàn)，顆粒的半徑服從均勻分布。 圖3顯示了花崗巖試樣在單軸壓縮試驗(yàn)和巴西劈裂試驗(yàn)下的應(yīng)力－應(yīng)變曲線，通過對(duì)比看出數(shù)值模擬得到的宏觀參數(shù)與實(shí)驗(yàn)參數(shù)基本一致。 校準(zhǔn)得到的巖石顆粒微觀參數(shù)如下:密度為2 500 kg/m3、半徑為0.9 ～1.5 mm、局部阻尼為0、孔隙度為0.16、有效模量為2.81 GPa、摩擦系數(shù)為0.7、粘結(jié)抗拉強(qiáng)度為85 MPa、粘結(jié)內(nèi)聚力為170 MPa、剛度比kn/ks為1.5、摩擦角φ 為30°。

圖3 微觀參數(shù)標(biāo)定曲線圖

此外，在記錄爆炸過程方面，由于目前無直接的方法得到炮孔壁單元速度、時(shí)程曲線，只能依托數(shù)值模擬再現(xiàn)爆炸過程。 而LS－DYNA 在解決非線性動(dòng)力問題和有效模擬爆炸過程方面具有突出的優(yōu)勢。因此，采用LS－DYNA 在爆炸試驗(yàn)基礎(chǔ)上建模，并根據(jù)實(shí)驗(yàn)構(gòu)建LS－DYNA 模型，其參數(shù)根據(jù)文獻(xiàn)[15]校準(zhǔn)。 爆炸應(yīng)力波導(dǎo)出及施加示意圖如圖4 所示。其中，炮孔模型及監(jiān)測點(diǎn)位置如圖4(a)所示。 爆炸產(chǎn)物壓力與體積的關(guān)系利用狀態(tài)方程JWL(Jones－Wilkins－Lee)描述，由式(3)表示為

圖4 爆炸應(yīng)力波導(dǎo)出及施加示意圖

式中P 為爆轟產(chǎn)物的壓力，GPa；V 為爆轟產(chǎn)物的相對(duì)體積；E0為單位體積的初始比內(nèi)能，GPa；A′、ω 為與炸藥類型有關(guān)的材料常數(shù)；B、R1、R2均為與爆炸物類型有關(guān)的材料常數(shù)。

炸藥材料參數(shù)如下:密度ρ ＝1 320 kg/m3、A′＝586 GPa、B ＝21.6 GPa、R1＝5.81、R2＝1.77、ω ＝0.282、E0＝7.38 GPa。

采用校準(zhǔn)過的微觀參數(shù)建立與現(xiàn)場試驗(yàn)相同的PFC2D 圓形巖石模型，刪除中間炮孔顆粒，PFC 模型中爆炸應(yīng)力波導(dǎo)出如圖4(a)所示；將單元速度均勻地施加在爆破孔邊界顆粒上，如圖4(b)所示。 圖5 顯示了LS－DYNA 得到的爆炸后爆破孔周圍單元的速度時(shí)程曲線。 具體軟件操作流程圖如圖6所示。

圖5 爆炸后爆破孔周圍單元的速度時(shí)程曲線圖

圖6 LS－DYNA/PFC 聯(lián)合方法流程圖

1.3 有效性驗(yàn)證

基于LS－DYNA/PFC2D 聯(lián)合分析的爆破過程中花崗巖裂紋擴(kuò)展模式如圖7(a) ～(e)所示，圖7(f)為文獻(xiàn)[6]中的電鏡實(shí)驗(yàn)裂紋圖。 對(duì)比選取的最終爆破與實(shí)驗(yàn)結(jié)果，可以看出兩者吻合度較高。在炸藥爆炸后，所引起的沖擊波會(huì)立即作用在爆破孔壁上，爆炸釋放的能量表現(xiàn)為從爆腔向外移動(dòng)的沖擊波荷載，炮孔周邊的巖石動(dòng)態(tài)抗壓強(qiáng)度遠(yuǎn)小于沖擊波強(qiáng)度，爆破孔周圍的圓形區(qū)域產(chǎn)生了致密的拉伸和剪切裂紋，爆破孔附近首先出現(xiàn)了破碎區(qū)，消耗了沖擊波大部分能量，破碎區(qū)內(nèi)巖石主要發(fā)生剪切破壞，在PFC2D 中顯示為紅色裂紋，如圖7(a)所示，此時(shí)炮孔周圍的裂紋數(shù)量增長迅速，沖擊波隨著與爆破孔距離的增加而衰減為應(yīng)力波，強(qiáng)度低于巖石的動(dòng)態(tài)抗壓強(qiáng)度，壓應(yīng)力波將炮孔周邊巖石向外推，破碎區(qū)外形成了拉伸區(qū)，巖石在環(huán)箍方向上產(chǎn)生拉伸應(yīng)力，遠(yuǎn)大于巖石的動(dòng)態(tài)抗拉強(qiáng)度，拉伸區(qū)巖石主要發(fā)生拉伸破壞，在PFC2D 中顯示為藍(lán)色裂紋。由于應(yīng)力波相對(duì)于沖擊波波速降低，且只在幾個(gè)方向產(chǎn)生徑向拉伸裂紋，裂紋數(shù)量增速放緩，此階段如圖7(b)所示，根據(jù)裂紋時(shí)間增長曲線可知此階段只發(fā)生拉伸破壞，當(dāng)沖擊波傳播至自由邊界處發(fā)生反射，形成拉伸波，并在巖石的邊界處產(chǎn)生剝落裂紋，拉伸裂紋數(shù)目迅速增加，隨后巖石周圍的剝落裂紋與炮孔周圍的拉伸徑向裂紋部分貫通，形成最后的裂紋分布。 由圖7(e)和(f)對(duì)比可以看出，LSDYNA/PFC 聯(lián)合模型的模擬結(jié)果與試驗(yàn)結(jié)果吻合良好，模型能夠準(zhǔn)確模擬巖石的損傷和破壞，因此該模擬方法可以有效地應(yīng)用于研究地應(yīng)力與爆破荷載聯(lián)合作用下巖石的裂紋擴(kuò)展規(guī)律。

圖7 巖石模擬及實(shí)驗(yàn)爆炸沖擊波下裂紋擴(kuò)展圖

2 地應(yīng)力對(duì)巖體爆破的影響分析

2.1 建立模型

通過1.3 節(jié)可以看出， LS－DYNA/PFC 聯(lián)合模型能夠準(zhǔn)地確模擬巖石在爆破荷載下的損傷破壞情形，通過建立方形巖體模型分析地應(yīng)力對(duì)爆破的影響，巖石的長、寬均為300 cm。 對(duì)模型施加應(yīng)力，考慮了7 種不同的分析工況，將x、y 方向的壓力Px和Py施加于模型的邊界，如圖8 所示，應(yīng)力施加工況見表1，同時(shí)考慮到巖石尺寸，采用2 倍直徑的炸藥進(jìn)行了研究。 通過銅管包圍，記錄LS－DYNA 炮孔周邊的單元速度，施加在方向模型上。 在建立的數(shù)值模型中創(chuàng)建6 個(gè)分區(qū)，觀察巖石在沖擊波及應(yīng)力波下的裂紋發(fā)育情況，詳細(xì)記錄裂紋數(shù)量的變化情況，如圖8(b)所示。 1、2 區(qū)分別為貼近炮孔區(qū)與水平和豎直方向成90°的扇形區(qū)域，其半徑為25 mm；緊接3、4 區(qū)位于1、2 區(qū)以外的扇形區(qū)域，其半徑為100 mm；而最遠(yuǎn)處為5、6 分區(qū)，其半徑為150 mm。

表1 不同條件下的原位應(yīng)力表

圖8 PFC 模型示意圖

2.2 單、雙軸應(yīng)力下爆炸分析

單軸應(yīng)力條件下爆破作用最終的微裂紋圖像如圖9 所示。 爆破荷載下徑向裂紋在水平方向上的擴(kuò)展長度隨著應(yīng)力的增大而減小，垂直于主應(yīng)力方向的裂紋長度隨著應(yīng)力的增大而增大。 圖10 顯示了巖石裂紋監(jiān)測的分區(qū)情況，以及各個(gè)分區(qū)裂紋數(shù)目在不同應(yīng)力下的對(duì)比情況。 在數(shù)值模型中創(chuàng)建6 個(gè)分區(qū)詳細(xì)記錄裂紋數(shù)量變化。 對(duì)比不同圍壓下各個(gè)分區(qū)裂紋數(shù)目，由于5、6 區(qū)未產(chǎn)生裂紋，此區(qū)域不做記錄，剪切裂紋只在靠近炮孔的1、2 區(qū)產(chǎn)生。 通過圖10(a)和(b)可知，圍壓對(duì)1、2 區(qū)影響較小，其裂紋數(shù)目高度重疊，由于炮孔近端沖擊應(yīng)力遠(yuǎn)大于地應(yīng)力，地應(yīng)力對(duì)破碎區(qū)域的影響較小。 如圖9(c)所示，與應(yīng)力垂直方向的4 區(qū)中裂紋數(shù)目隨著應(yīng)力增加而減小，表現(xiàn)為應(yīng)力對(duì)裂紋的產(chǎn)生有抑制作用，根據(jù)圖10(d)在不同應(yīng)力下裂紋數(shù)目對(duì)比可知，平行于應(yīng)力方向的裂紋在一定范圍內(nèi)會(huì)增加，應(yīng)力對(duì)裂紋的產(chǎn)生有促進(jìn)作用。 距離炮孔較遠(yuǎn)的3、4 區(qū)域裂紋顯然受應(yīng)力影響，地應(yīng)力對(duì)裂縫數(shù)目的影響隨著距爆破孔距離的增加而逐漸增加。

圖9 單軸應(yīng)力下微裂紋擴(kuò)展分布圖

圖10 不同區(qū)域單軸應(yīng)力下微裂紋數(shù)量對(duì)比圖

圖11 為雙軸應(yīng)力條件下爆炸作用的裂縫最終分布情況。 隨著應(yīng)力的增加裂，紋長度逐漸變短，圍壓對(duì)拉伸裂紋影響較明顯。 由于雙軸應(yīng)力下裂紋分布對(duì)稱，所以不再分區(qū)監(jiān)測，而采用總裂紋數(shù)量進(jìn)行對(duì)比，如圖12 所示，微裂紋數(shù)量隨著應(yīng)力的增加明顯下降。 從裂紋數(shù)目及裂紋形狀可以看出，原位應(yīng)力對(duì)裂紋的發(fā)展有抑制作用。

圖11 雙軸應(yīng)力下微裂紋最終分布圖

圖12 雙軸應(yīng)力下炮孔微裂紋數(shù)量演變圖

雙軸應(yīng)力下破碎區(qū)的形狀是圓形的，炮孔周圍裂紋主要由剪切裂紋和拉伸裂紋組成，隨著距離的增加，幾個(gè)主要裂紋從破碎區(qū)周邊徑向延伸，各個(gè)方向裂紋長度幾乎相同，隨著應(yīng)力的增加，最遠(yuǎn)端徑向裂紋明顯長度減小，裂紋數(shù)量減少明顯。 根據(jù)圖12中4 種應(yīng)力下裂紋數(shù)量演變過程對(duì)比可知，爆炸前期裂紋數(shù)量增長趨勢相差不大，即原位應(yīng)力對(duì)爆破孔周圍的破碎區(qū)影響不大，隨著距離炮孔距離的增加，應(yīng)力對(duì)裂紋數(shù)量的影響越來越大，遠(yuǎn)處徑向裂紋的延伸受到強(qiáng)烈抑制。

為了研究在等雙軸應(yīng)力和爆破荷載作用下巖石中的應(yīng)力分布，在炮孔的右側(cè)40 mm 處設(shè)置測量圓，詳細(xì)分析測量圓中的動(dòng)態(tài)應(yīng)力響應(yīng)。 巖石中某點(diǎn)徑向應(yīng)力和切向應(yīng)力計(jì)算公式分別由式(4)和(5)表示為

式中σn和τn分別為正應(yīng)力和切應(yīng)力，MPa；σx和σy分別為x、y 軸方向的應(yīng)力分量，MPa；τxy＝τyx為切向應(yīng)力分量， MPa； n 為平面AB 的法線；l 和m 為應(yīng)力主軸的兩個(gè)方向余弦； θ 為與y 軸正向的夾角，°。任一點(diǎn)徑向應(yīng)力σr和切向應(yīng)力σθ均可由式(4) 和(5) 計(jì)算得出。

在PFC2D 中，測量圓可以記錄應(yīng)力分量的時(shí)間曲線(x、y、τ)，測量圓中記錄的應(yīng)力包括靜態(tài)應(yīng)力和動(dòng)態(tài)應(yīng)力，各應(yīng)力分量分布示意如圖13 所示。PFC 程序中規(guī)定測量圓應(yīng)力值為負(fù)時(shí)表示巖體處于壓縮狀態(tài)，反之處于拉伸狀態(tài)；測量圓的應(yīng)變速率是根據(jù)顆粒間的距離變化來測量的，應(yīng)變速率為負(fù)時(shí)表示沿受力方向巖體內(nèi)顆粒的平均距離變小，反之表示顆粒間平均距離變大。 通常監(jiān)測點(diǎn)一般會(huì)出現(xiàn)兩個(gè)峰值，即壓縮應(yīng)力峰值與之后的拉伸應(yīng)力峰值，峰值強(qiáng)度由應(yīng)力波的強(qiáng)度與初始靜態(tài)應(yīng)力決定，圖14(a)和(b)分別顯示了在等雙軸應(yīng)力為0、5、10 和15 MPa 時(shí)徑向應(yīng)力－時(shí)間和切向應(yīng)力－時(shí)間曲線。從圖11 可以看出，在等雙軸原位應(yīng)力和爆破載荷耦合下，徑向壓應(yīng)力隨原位應(yīng)力的增加而上升，相反環(huán)箍向張拉應(yīng)力隨原位應(yīng)力的增加而減小，過高的初始應(yīng)力導(dǎo)致拉伸應(yīng)力峰值低于巖石動(dòng)態(tài)抗拉強(qiáng)度，不會(huì)造成監(jiān)測區(qū)域的張拉破壞。 由此得出，在等雙軸應(yīng)力下拉伸裂紋被抑制，并且隨著原位應(yīng)力的升高，抑制效果也隨之增強(qiáng)。

圖13 巖石平面應(yīng)力分布示意圖

圖14 不同地應(yīng)力下的測點(diǎn)應(yīng)力－時(shí)間曲線圖

3 結(jié)論

文章采用LS－DYNA/PFC 相結(jié)合的方法研究爆炸沖擊波下裂紋變化，模擬了沖擊波引起的破碎區(qū)及裂紋區(qū)，分析了地應(yīng)力在巖石爆破過程中的影響特性，得到如下主要結(jié)論:

(1) 在雙軸等圍壓情況，微裂紋數(shù)目曲線在爆炸初期階段高度重合，圍壓對(duì)破碎區(qū)影響較小，隨著距離增加，裂紋數(shù)量差異將越來越明顯，圍壓對(duì)拉伸區(qū)裂紋的擴(kuò)展影響也較大；炮孔附近的徑向、環(huán)向壓應(yīng)力均隨地應(yīng)力的升高而升高；預(yù)壓應(yīng)力增強(qiáng)了壓縮效果，削弱了炮孔徑向方向的張力效應(yīng)，阻礙了裂紋擴(kuò)展。

(2) 在雙軸不等圍壓情況下，主應(yīng)力的方向和大小并不會(huì)影響炮孔近端破碎區(qū)裂紋數(shù)目變化，僅影響較遠(yuǎn)端裂紋帶的發(fā)展。 裂紋發(fā)展方向與單軸壓力方向平行時(shí)，單軸壓力促進(jìn)了裂紋的形成；裂紋方向與單軸壓力方向垂直時(shí)，單軸壓力抑制了裂紋的形成，且隨著壓力的增加，其效果越明顯。