周文海，梁 瑞，陳金林，朱 冕，陳鵬輝，樓曉明，王敦繁

（1. 蘭州理工大學(xué)石油化工學(xué)院，甘肅 蘭州 730050；2. 紫金礦業(yè)集團(tuán)股份有限公司紫金山金銅礦，福建 龍巖 364200；3. 福州大學(xué)爆炸技術(shù)研究所，福建 福州 350116；4. 浙江大學(xué)海洋學(xué)院，浙江 舟山 316021）

選取合理的臺階逐孔微差爆破延期時(shí)間，可以降低地表質(zhì)點(diǎn)振動，改善爆破質(zhì)量，節(jié)省工程投入成本；自20 世紀(jì)中葉，通過控制微差時(shí)間降低爆破振動危害和改善爆破質(zhì)量方面有了大量科學(xué)研究[1]。田振農(nóng)等[2]采用時(shí)-頻分析方法，對隧道爆微差破振動的一般特征進(jìn)行了研究，發(fā)現(xiàn)：由于雷管段數(shù)的限制，起爆網(wǎng)絡(luò)經(jīng)常被分成幾個(gè)時(shí)段，造成振動波無規(guī)律疊加，以至出現(xiàn)多個(gè)峰值現(xiàn)象；但是如果技術(shù)措施合理，能使爆破振動速度峰值顯著降低，在此基礎(chǔ)上根據(jù)雷管起爆延時(shí)精度高的特點(diǎn)并借鑒干擾減振的思路，提出的錯相減振機(jī)理，對微差起爆振動效應(yīng)起到了有效控制。Johansson 等[3]通過微差控制小炮實(shí)驗(yàn)對孔間沖擊波相互作用進(jìn)行研究，通過大量研究數(shù)據(jù)表明，微差時(shí)間選?。?～1.1 W ms）時(shí)降振效果明顯優(yōu)于2 W ms 及以上，其中W 為最小抵抗線（單位為m）。龔敏等[4]以南方某城市隧道工程為背景，進(jìn)行現(xiàn)場單孔單自由面爆破實(shí)驗(yàn)，獲得不同藥量-微差時(shí)間振動曲線，再利用MATLAB 程序，根據(jù)每段不同延時(shí)范圍，將兩單孔曲線按相鄰段起爆的多個(gè)微差間隔進(jìn)行不同振動疊加，選擇其中最大振速的合成曲線與單孔振動曲線按下一相鄰孔的微差間隔進(jìn)行新的振動合成，最終得到微差掏槽爆破后的累積疊加曲線，以此與現(xiàn)場實(shí)測振動曲線比較，得出微差起爆間隔較大的合成振速并不一定比小間隔起爆的振速小，合成振速是否超標(biāo)取決于起爆藥量(或相應(yīng)振動曲線)、微差間隔時(shí)間、規(guī)定振速值三者的量化關(guān)系。Anderson 等[5]運(yùn)用統(tǒng)計(jì)分析理論，對大量微差爆破工程實(shí)例資料進(jìn)行分析研究，最終利用灰度圖反演了微差時(shí)間與振動波主頻之間關(guān)系，以此劃分出了不同微差時(shí)間對于地表振動的影響程度和范圍。楊年華等[6]將現(xiàn)場試爆單孔藥包引起的振動波作為基波，按照不同比藥量的比例系數(shù)對基波進(jìn)行折算得到實(shí)際施工振動波譜，然后按不同微差時(shí)間利用振動波線性疊加理論對折算后的振動波譜合成，以此達(dá)到對峰值振動強(qiáng)度預(yù)測和對最佳微差時(shí)間進(jìn)行判定的目的。本文中，先通過有限元法確定自然狀態(tài)下二維邊坡潛在滑動面以及靜態(tài)安全系數(shù)，再基于二維靜力模型重新建立三維爆破模型，反演沖擊載荷作用下潛在滑動面受力情況，并提取相應(yīng)的滑面單元應(yīng)力數(shù)值，結(jié)合傳統(tǒng)極限平衡法繪制出沖擊載荷作用下邊坡穩(wěn)定性系數(shù)曲線，最終確定出最佳孔間降振微差時(shí)間。

1 沖擊載荷作用下安全系數(shù)計(jì)算方法

通常有兩種方法用于邊坡穩(wěn)定性分析，一種為極限平衡法，另一種則是有限元法。極限平衡法將滑坡體看作理想塑性材料的均質(zhì)剛性體，計(jì)算過程對邊界條件大大簡化，且完全不考慮巖土體本身應(yīng)力-應(yīng)變關(guān)系，因而很難真實(shí)反映邊坡失穩(wěn)狀態(tài)下位移場和應(yīng)力場的變化情況；而有限元法則考慮了巖土體應(yīng)力-應(yīng)變關(guān)系，且不需要預(yù)先設(shè)定滑動面，能夠真實(shí)有效地反映巖土體邊坡失穩(wěn)狀態(tài)。因此，本文中先利用有限元法確定潛在滑動面，提取滑動面單元上的應(yīng)力 σx、 σy、 τxy，再結(jié)合極限平衡法求 得邊坡體在沖擊載荷下的動力穩(wěn)定性系數(shù)，如圖1所示。

極限平衡法以摩爾-庫倫抗剪強(qiáng)度理論為基礎(chǔ)，對潛在滑動面進(jìn)行條塊劃分，在各條塊上建立平衡方程，將潛在滑動面上抗滑力和滑動力的比值f 與定義為安全系數(shù)：

式中： τf為巖土體抗剪強(qiáng)度； τ為巖土體剪應(yīng)力。

運(yùn)用有限元法確定潛在滑動面時(shí)需要定義巖土體本構(gòu)模型，而ANSYS 程序中自帶多種屈服破壞準(zhǔn)則，巖土體材料通常選用摩爾-庫倫準(zhǔn)則或廣義米塞斯準(zhǔn)則（D-P 準(zhǔn)則）[7]。在考慮靜水壓力情況下，將米塞斯準(zhǔn)則進(jìn)行推廣就轉(zhuǎn)換成D-P 準(zhǔn)則：

式中： σ1、 σ2、 σ3分別為第一、二、三主應(yīng)力； I1為應(yīng)力張量第一不變量； J2為應(yīng)力偏張量第二不變量； α、k為模型材料與內(nèi)摩擦角 φ′和凝聚力 c′相關(guān)系數(shù)。

基于有限元計(jì)算結(jié)果，可提取滑面單元上的應(yīng)力 σx、 σy、 τxy。 切向應(yīng)力τ 和法向應(yīng)力 σn為：

式中： α 為 x軸與所求滑面單元切線方向夾角。

由摩爾-庫倫準(zhǔn)則，可得所求滑面單元抗剪強(qiáng)度為：

式中： c′為巖土體凝聚力； φ′為巖土體內(nèi)摩擦角。

將整個(gè)邊坡潛在滑動面單元按式（1）進(jìn)行積分運(yùn)算，可得邊坡整體安全系數(shù)公式為：

式中： l為潛在滑坡面弧段長度。

2 按實(shí)體模型進(jìn)行數(shù)值模擬

2.1 按實(shí)體模型建模

以西北某露天采場臨時(shí)邊坡為原型進(jìn)行數(shù)值模擬，建模所用材料參數(shù)以爆破施工現(xiàn)場的具體情況為準(zhǔn)。其中，模擬過程所采用的巖石力學(xué)性能源于該礦區(qū)地質(zhì)資料中的爆區(qū)附近采樣點(diǎn)數(shù)據(jù)，爆破工藝參數(shù)和模型尺寸與施工情況一致，其邊坡高度24 m，坡面角70°，孔徑150 mm，孔深14 m，堵塞長度5 m，底盤抵抗線根據(jù)鉆孔作業(yè)安全條件 W=H cotα+B以臺階高度12 m 來設(shè)計(jì)，取W=7 m（H 為臺階高度、 α為臺階坡面角；B 為鉆孔中心至坡頂線安全距離取2.5 m），孔排間距選用5 m×6 m 的孔網(wǎng)參數(shù)。巖石材料選取典型的彈塑性材料，即采用MAT_PLASTIC_KINEMATIC 模型；堵塞物選用土壤材料MAT_SOIL_AND_FOAM 模型；炸藥選取2#巖石乳化炸藥，其材料類型為HIGH_EXPLOSIVE_BURN，炸藥狀態(tài)方程選取不考慮炸藥產(chǎn)物成分的EOS_JWL 方程。為提高模擬運(yùn)算效率，在不影響計(jì)算精度和分析結(jié)果的基礎(chǔ)上：首先通過ANSYS 軟件建立自然狀態(tài)下二維實(shí)體邊坡靜力模型，運(yùn)用有限元折減法確定潛在滑動面；然后依據(jù)二維靜力模型，重新建立同性質(zhì)同尺寸的三維爆破模型，通過LS-DYNA 軟件進(jìn)行逐孔微差爆破動力分析；整個(gè)模擬過程所用爆破參數(shù)、巖土力學(xué)參數(shù)和邊坡結(jié)構(gòu)參數(shù)按照現(xiàn)場實(shí)際取值。模擬過程炸藥爆炸時(shí)的JWL 狀態(tài)方程[8]為：

式中：A、B、 R1、 R2、 ω為 材料常數(shù)，p 為壓力，V 為相對體積， E為初始比內(nèi)能。

巖石主要力學(xué)參數(shù)分別為：巖土密度為2.78 g/cm3，彈性模量E=18.23 GPa，抗壓強(qiáng)度σ0=76.2 MPa，泊松比μ=0.23，內(nèi)聚力為17.6 MPa，內(nèi)摩擦角φ=33°36'。炸藥主要參數(shù)分別為：密度為0.95 g/cm3，爆速為3.6 km/s，爆轟壓力為3.61 GPa，A=47.6 GPa，B=529 MPa，R1=3.5，R2=0.9，E=4.5 GPa，V0=1 cm3。

鄭穎人等[7]、Griffiths 等[9]通過對比研究指出，有限元模擬時(shí)模型邊界范圍取值對于計(jì)算結(jié)果的影響較極限平衡法更敏感，且給出較理想的模型范圍取值為：邊坡上下邊界高度應(yīng)超出坡高2 倍以上，右邊界距坡頂線距離取坡高2.5 倍左右長度，左邊界距坡腳處距離取坡高1.5 倍左右長度；因此，將三維實(shí)體邊坡模型進(jìn)行拓展，上下邊界取50 m，左邊界距坡腳距離取36 m，右邊界距坡頂線距離取60 m，邊坡厚度取50 m，模型三維示意圖和參數(shù)平面圖如圖2～3 所示。

圖 2 三維有限元計(jì)算模型Fig. 2 3D finite element model of slope

圖 3 模型參數(shù)Fig. 3 Model parameters

2.2 單元類型選取、網(wǎng)格劃分及施加載荷

二維建模采用四節(jié)點(diǎn)的PLANE82 單元和理想彈塑性模型（D-P 模型）；三維建模靜力分析時(shí)采用八節(jié)點(diǎn)隱式實(shí)體單元SOLID 185，動力分析時(shí)采用顯示分析和八節(jié)點(diǎn)SOLID 164 實(shí)體單元，材料定義為理想彈塑性模型。模型網(wǎng)格劃分嚴(yán)格按照二維狀態(tài)四邊形、三維狀態(tài)六面體進(jìn)行，二維、三維模型劃分單位數(shù)分別為4 521 和170 958。在施加約束和載荷方面：二維模型邊坡面和頂部定義為自由邊界，左右邊界面施加位移約束，運(yùn)算過程一直保持重力加載；三維模型邊坡面和頂部同樣定義為自由邊界，其余均為透射邊界，運(yùn)算過程保持動力松弛加載預(yù)應(yīng)力且全程施加重力載荷。二維、三維網(wǎng)格劃分如圖4～5 所示。

圖 4 二維模型網(wǎng)格劃分Fig. 4 Mesh generation of 2D

圖 5 三維模型網(wǎng)格劃分Fig. 5 Mesh generation of 3D

2.3 自然狀態(tài)下靜態(tài)安全系數(shù)的確定

采用強(qiáng)度折減法判定邊坡是否失穩(wěn)的充要條件為：滑動面貫通于塑性應(yīng)變區(qū)、滑動面上節(jié)點(diǎn)位移或塑性應(yīng)變發(fā)生突變、有限元方程組無解且計(jì)算結(jié)果不收斂[10-15]。依據(jù)有限元強(qiáng)度折減原理，先選取初始折減系數(shù)為1，將巖土體強(qiáng)度參數(shù)（內(nèi)摩擦角、黏聚力）進(jìn)行折減，當(dāng)?shù)? 次折減運(yùn)算時(shí)（ f=2.750），潛在滑動面處于臨界破壞狀態(tài)，直到安全系數(shù)取 f=2.780 時(shí)滑動面上節(jié)點(diǎn)位移發(fā)生突變且貫通于塑性應(yīng)變區(qū)，有限元計(jì)算結(jié)果不收斂。具體塑性應(yīng)變?nèi)鐖D6 所示，折減過程參數(shù)變化見表1。

將強(qiáng)度折減法計(jì)算得到的潛在滑動面按照傳統(tǒng)極限平衡法進(jìn)行條塊劃分，并且計(jì)算條塊滑弧每個(gè)單元長度和滑弧面上的切線與x 軸的夾角θ，進(jìn)行邊坡穩(wěn)定性驗(yàn)證分析。通過計(jì)算分析發(fā)現(xiàn)：二維有限元臨界安全系數(shù)為2.750，三維有限元臨界安全系數(shù)為2.610；傳統(tǒng)極限平衡法中瑞典條分法確定的臨界安全系數(shù)為2.525，簡化Bishop 法確定的臨界安全系數(shù)為2.505，Jonbu 確定的臨界安全系數(shù)為2.679。因此，有限元法確定的臨界安全系數(shù)與傳統(tǒng)極限平衡法確定的臨界安全系數(shù)相對誤差在8%以內(nèi)，計(jì)算結(jié)果較接近。各計(jì)算方式得到的安全系數(shù)見表2。

圖 6 塑性區(qū)域分布Fig. 6 Plastic regional distribution

表 1 折減算法Table 1 Strength reduction method

表 2 各種計(jì)算方式的安全系數(shù)Table 2 Safety factor for each calculation method

2.4 沖擊載荷作用下動態(tài)安全系數(shù)的確定

根據(jù)二維模型得到的邊坡滑動面，依據(jù)同等尺寸重新建立同等條件下的三維逐孔微差爆破模型，將邊坡潛在滑坡弧面進(jìn)行等長條分，取三維臺階中間剖面進(jìn)行網(wǎng)格細(xì)分，共等分為37 個(gè)單元；臺階面同排設(shè)置3 個(gè)炮孔微差起爆，孔間微差時(shí)間分別取0、17、25、42 和65 ms，如圖1～2 所示。利用ANSYS 軟件進(jìn)行顯式-隱式轉(zhuǎn)換，通過動力松弛分析，得到邊坡在重力作用下的原始應(yīng)力場；將ANSYS 隱式求解結(jié)果導(dǎo)出的動力松弛文件導(dǎo)入LS-DYNA 中進(jìn)行后續(xù)動力分析，3 個(gè)炮孔同時(shí)起爆某個(gè)時(shí)刻三維應(yīng)力云分布如圖7 所示。

圖 7 應(yīng)力分布Fig. 7 Stress distribution

對3 個(gè)炮孔同時(shí)起爆過程中模型質(zhì)點(diǎn)的應(yīng)力、位移進(jìn)行分析，發(fā)現(xiàn)：起爆初期炮孔壁壓力急劇增大，最大值可達(dá)1.36 GPa，直到約60 ms 孔壁壓力基本趨于穩(wěn)定，孔壁單元應(yīng)力明顯高于坡面應(yīng)力；坡腳處出現(xiàn)應(yīng)力集中現(xiàn)象，應(yīng)力明顯高于其他質(zhì)點(diǎn)，但應(yīng)力集中效應(yīng)瞬間消失后趨于穩(wěn)定。坡頂質(zhì)點(diǎn)位移波動幅值高于坡面和坡腳位移約3 倍，坡頂位移波動幅值分布于?0.35～0.25 cm，這主要是應(yīng)力波在自由面處的反射拉伸作用和尖端放大效應(yīng)所致。

圖 8 微差起爆合速度衰減規(guī)律Fig. 8 Resultant velocity attenuation regularity of millisecond detonating

依據(jù)樓曉明等[16]提出的逐孔微差爆破振動波峰值合速度-位移分布特征理論，可將爆破振動過程中地表質(zhì)點(diǎn)三軸振動峰值合速度， Vs,max作為地表質(zhì)點(diǎn)振動大小判據(jù)；因而在滑坡弧中間剖面提取其中某個(gè)單元觀測其合速度 Vs變化規(guī)律，通過觀察發(fā)現(xiàn)17、25、42 和65 ms 不同孔間微差延期起爆時(shí)對應(yīng)的 Vs,max分別為9.08、14.3、6.93 和9.12 cm/s，說明選用42 ms 孔間微差延期起爆時(shí)間能夠起到更好的降振效果。4 種微差起爆 Vs衰減規(guī)律如圖8 所示。

圖 9 邊坡穩(wěn)定性系數(shù)曲線Fig. 9 Slope stability coefficient curve

2.5 安全系數(shù)曲線確定微差時(shí)間

利用后處理軟件提取3 個(gè)炮孔同時(shí)起爆時(shí)滑面單元爆炸過程中的應(yīng)力 σx、 σy、 τxy；在1～300 ms 運(yùn)算時(shí)間段內(nèi)，每隔3 ms 輸出一次，共計(jì)100 組數(shù)據(jù)，利用式（5）～（8）求解動態(tài)安全系數(shù)，并繪制出邊坡穩(wěn)定性系數(shù)曲線，如圖9 所示。

由圖9 可知，藥包起爆后在0～21 ms 時(shí)間段，內(nèi)應(yīng)力波還未集中作用于邊坡潛在滑動面，只有少部分應(yīng)力波抵達(dá)該部位，所以該時(shí)間段內(nèi)邊坡動態(tài)安全系數(shù)基本和自然狀態(tài)下靜態(tài)穩(wěn)定性系數(shù)保持一致。隨著時(shí)間推移，在爆炸沖擊載荷作用下邊坡動態(tài)穩(wěn)定性發(fā)生較大波動，當(dāng)t=48 ms 時(shí)動態(tài)穩(wěn)定性數(shù)達(dá)到最大值 f =7.623，當(dāng)t=102 ms 時(shí)動態(tài)穩(wěn)定性系數(shù)達(dá)到最小值 f=1.375，直到240 ms 左右趨于穩(wěn)定。

從先爆炮孔產(chǎn)生的沖擊載荷集中作用于邊坡潛在滑動面開始，到該潛在滑動面穩(wěn)定性系數(shù)達(dá)到最大值（最穩(wěn)定狀態(tài)）為止，將該時(shí)間間隔步長作為后續(xù)炮孔起爆的孔間最佳降振微差時(shí)間較合理。通過本文選取的實(shí)際邊坡爆破模型進(jìn)行孔間同時(shí)起爆數(shù)值模擬發(fā)現(xiàn)，21 ms 時(shí)沖擊載荷集中作用于邊坡潛在滑動面，48 ms 時(shí)動態(tài)穩(wěn)定性系數(shù)達(dá)到最大值；因此，對于該模型選用孔間微差爆破技術(shù)時(shí)，最佳孔間降振微差時(shí)間應(yīng)為48 ms。

3 測振實(shí)驗(yàn)選取降振微差時(shí)間

在相同露天采場臨時(shí)邊坡平臺，進(jìn)行了4 次與模擬情況相似的等比例不同段別雷管組合的小炮爆破測振實(shí)驗(yàn)，4 次實(shí)驗(yàn)所選用的部分爆破工藝參數(shù)取值與模擬參數(shù)比約為1∶2。小炮實(shí)驗(yàn)的孔徑80 mm，孔深7 m，堵塞2.5 m，底盤抵抗線W=4 m，孔排間距選用3 m×4 m 的孔網(wǎng)參數(shù)，裝藥結(jié)構(gòu)為無間隔耦合裝藥。4 次實(shí)驗(yàn)排間統(tǒng)一選用3 段25 ms 延期雷管，孔間分別采用2、3、4、5 段雷管控制微差時(shí)間，即組成排、孔間（25 ms，17 ms）、（25 ms，25 ms）、（25 ms，42 ms）、（25 ms，65 ms）的微差起爆方式；試爆時(shí)，采用爆破測振儀對距爆區(qū)空間距離約270 m 范圍進(jìn)行測振，每次測振區(qū)布置4 個(gè)臨近監(jiān)測點(diǎn)，分別記為A、B、C、D，測振所得峰值速度數(shù)據(jù)和4 次監(jiān)測中A 測點(diǎn)合速度衰減規(guī)律見表3和如圖10 所示。

通過對4 次微差試爆區(qū)16 個(gè)測點(diǎn)測得的三軸峰值合速度數(shù)值分析發(fā)現(xiàn)：孔間選取17 ms 微差起爆時(shí)，4 個(gè)監(jiān)測點(diǎn)測得的峰值合速度最大值為1.294 1 cm/s，最小值為1.067 7 cm/s；孔間選取25 ms 微差起爆時(shí)，4 個(gè)監(jiān)測點(diǎn)測得的峰值合速度最大值為1.702 1 cm/s，最小值為1.387 8 cm/s；孔間選取42 ms微差起爆時(shí)，4 個(gè)監(jiān)測點(diǎn)測得的峰值合速度最大值為0.630 6 cm/s，最小值為0.503 1 cm/s；孔間選取65 ms 微差起爆時(shí)，4 個(gè)監(jiān)測點(diǎn)測得的峰值合速度最大值為0.943 7 cm/s，最小值為0.865 5 cm/s。測試結(jié)果顯示，孔間選取42 ms 微差時(shí)間延期起爆時(shí)測點(diǎn)三軸峰值合速度幅值明顯低于其他3 種情況，依據(jù)峰值合速度-位移分布特征理論，可證實(shí)4 種測振實(shí)驗(yàn)中孔間微差時(shí)間取42 ms 時(shí)引起的地表質(zhì)點(diǎn)振動最小。

表 3 不同孔間微差時(shí)間下監(jiān)測點(diǎn)峰值速度Table 3 Monitoring point’s peak speed at different hole millisecond time

圖 10 不同孔間微差起爆時(shí)測點(diǎn)A 合速度衰減規(guī)律Fig. 10 Resultant velocity attenuation regularity of measuring point A with different millisecond detonating

4 結(jié) 論

(1)研究爆破施工對邊坡穩(wěn)定性影響時(shí)，可先運(yùn)用有限元法建立二維靜態(tài)邊坡模型確定潛在滑動面，然后在已確定的潛在滑動面基礎(chǔ)上建立三維爆破模型，同時(shí)采用LS-DYNA 軟件對沖擊載荷作用下邊坡穩(wěn)定性進(jìn)行動力分析。

(2)在模擬結(jié)果中，提取沖擊載荷作用下邊坡滑動面單元應(yīng)力，將該應(yīng)力代入傳統(tǒng)極限平衡法計(jì)算公式，繪制出邊坡穩(wěn)定性系數(shù)曲線；通過對該曲線分析發(fā)現(xiàn)，可將沖擊載荷集中作用于邊坡滑動面開始到該滑動面穩(wěn)定性系數(shù)達(dá)到最大值之間的時(shí)間步長作為孔間最佳降振微差時(shí)間。

(3)通過對某個(gè)實(shí)際邊坡穩(wěn)定性系數(shù)曲線理論分析發(fā)現(xiàn)，最佳降振微差時(shí)間為48 ms，而通過不同微差時(shí)間控制的數(shù)值模擬和現(xiàn)場測振實(shí)驗(yàn)得到最佳降振微差時(shí)間為42 ms；該理論分析結(jié)果與模擬和實(shí)驗(yàn)結(jié)果較接近，說明采用邊坡動態(tài)穩(wěn)定性系數(shù)曲線確定的最佳降振微差時(shí)間較可靠，可為類似工程項(xiàng)目提供理論依據(jù)。