爆炸載荷下順層臺階邊坡漸進破壞規律數值分析1)

姚作強 * 馮 春 ** 劉天蘋

*(北京市工程咨詢有限公司，北京 100124)

?(中國科學院力學研究所流固耦合系統力學重點實驗室，北京 100190)

**(中國科學院大學工程科學學院，北京 100049)

臺階爆破在露天礦開采，水電邊坡、公路（鐵路）沿線邊坡開挖等方面廣泛應用。在臺階爆破過程中，爆破產生的應力波將對后方的高陡邊坡產生影響。近距離或大當量爆破過程中，在應力波反復作用下，邊坡內部的薄弱結構面將逐漸出現失穩貫通，最終可能導致邊坡出現整體性滑動。因此，臺階邊坡的爆破穩定性問題一直受到工程界及學術界的持續關注[1-7]。

數值模擬因其能揭示爆破作用下邊坡內部應力波的傳播規律及宏觀裂縫的漸進擴展規律，因此被廣泛應用于邊坡爆破穩定性的分析中。占紹祥等[8]和張卉等[9]重點分析了爆破開挖過程中Hoek–Brown 擾動系數與邊坡穩定性的對應關系，并指出采用擾動系數漸進弱化的方式，可較好反映爆破損傷對邊坡穩定的影響。王建國等[10]利用FLAC3D 軟件，通過施加實測地震波的方式，探討了布沼壩露天煤層西幫邊坡的動態規律。周后友等[11]借助ANSYS 分析了別礦露天邊坡在靜力和動力下的穩定性，給出了該邊坡在爆破載荷下的安全振速和臨界振速。吳燕開等[12]采用ABAQUS分析了爆炸載荷下巖石的振動效應，給出了巖質高邊坡在動力作用下的位移演化規律及振動規律。馬沖等[13]借助FLAC3D 探討了爆破載荷對邊坡振速規律的影響，并建立了爆破振速安全閾值。

目前，邊坡爆破穩定性的常用數值模擬方法包括有限元法、離散元法和有限差分法等。有限單元法通過在表征元上引入宏觀動力學本構，實現了爆炸載荷下邊坡彈塑性變形過程的精確模擬，但卻無法準確刻畫邊坡內部結構面的漸進失穩過程及邊坡整體失穩后的運動成災過程。離散元方法適用于模擬真實顆粒體系的碰撞運動過程，但在模擬連續介質彈塑性變形方面精度較差。連續–非連續單元方法(continuum–discontinuum element method，CDEM)[14-20]通過將連續介質數值模擬方法（如有限元法）及非連續介質數值模擬方法（如離散元法）進行有機結合，實現了地質體彈性、塑性、損傷、斷裂、破碎、運動、碰撞、堆積的全過程模擬。

總體而言，目前邊坡爆破穩定性的數值分析大都基于場地實測獲得的加速度（速度）時程曲線展開，沒有實現炸藥起爆、應力波傳播、邊坡振動破壞的全過程耦合。事實上，由于受場地條件及地層特性的影響，地表測得的振動時程曲線與通過巖體內部傳遞至坡體的應力波并不完全等價，因此上述將爆破過程與邊坡失穩過程割裂開來的分析方法，無法準確反映爆破對邊坡的損傷破壞作用。此外，目前的數值模擬大都以連續介質方法（如有限元法、有限差分法）為主，無法準確刻畫爆炸載荷下結構面巖體的漸進破裂過程及失穩后的運動成災過程。因此，本文將以CDEM的數值計算代碼為基礎，通過引入考慮應變率效應的巖體及結構面動力損傷破裂力學模型及爆生氣體的膨脹模型，編制C++ 數值計算代碼，直接模擬三維情況下爆破對順層巖質邊坡的影響，并重點探討爆破參數及結構面參數對邊坡穩定性的影響規律。

1 計算方法及本構模型

1.1 CDEM 計算方法

CDEM 是一種以廣義拉格朗日方程為基本框架的顯式動力學數值方法。該方法將基于連續與非連續介質的兩類數值方法進行深度融合，在能量層面實現了有限元法、離散元法及無網格法的統一，可模擬地質體及人工材料在靜、動力載荷下的彈性、塑性、損傷及破裂過程，以及破碎后散體的運動、碰撞、流動及堆積過程。CDEM 中每個單元均具有獨立的節點及獨立的面，相鄰兩個單元依靠粘接界面進行連接；利用單元表征材料宏觀連續介質特性，利用粘接界面表征材料的斷裂滑移等非連續介質特性（圖1）。

圖1 CDEM 中的斷裂描述Fig.1 Description of fracture in CDEM

CDEM 采用增量顯式迭代法進行動力問題的求解，主要迭代過程包括節點合力的求解及節點運動過程的求解兩個方面。

節點合力計算公式為

式中，F為節點合力，FE為節點外力，Fe為有限元單元變形貢獻的節點力，Fc為界面接觸在節點上貢獻的接觸力，Fr為瑞利阻尼在節點上貢獻的節點阻尼力，Fl為節點局部阻尼力。

節點運動計算公式為

式中，a為節點加速度，v為節點速度， ?u為節點位移增量，u為節點位移，m為節點質量，?t為計算時步。

基于式（1）及式（2）的交替計算即為顯式求解過程。

1.2 邊坡巖體本構模型

為了精確描述爆炸載荷下巖體的動力學行為及漸進破壞過程，本文提出了考慮應變率效應及損傷破裂效應的單元–界面雙重本構，采用考慮動力增強效應的Mohr–Coulomb 應變軟化模型對邊坡的實體單元進行描述，采用考慮應變率效應的斷裂能本構對巖體結構面進行描述。

（1）實體單元本構[21]

首先判斷當前應力狀態是否已經達到或超過Mohr–Coulomb 準則及最大拉應力準則，為

式中，fs為Mohr–Coulomb 屈服函數；ft為最大拉應力屈服函數；h為Mohr–Coulomb 屈服函數與最大拉應力屈服函數的角平分線（用于判定失效模式）；σ1和σ3為最小主應力及最大主應力；σt(t)，c(t)，?為當前步抗拉強度、黏聚力及內摩擦角；αp，σp，N?為常數；η 及ξ 分別為拉伸應變率增強因子及剪切應變率增強因子。

η 及ξ 可表述為

式中，θ˙ 及γ˙ 為體應變率及等效剪應變率，e1及e2為單元應變率效應指數，本文取e1=e2=1/3。˙

θ及 ˙γ可表述為

式（3）中，如果fs≥0 且h≤0，單元發生剪切塑性變形；如果ft≥0 且h>0，則發生拉伸塑性變形。

單元發生塑性變形后，采用式（6）對黏聚力及抗拉強度進行折減，為

式中，Δt為計算步長，σt0和c0為抗拉強度及黏聚力的初始值，εp和εlim為當前時刻拉伸塑性應變及拉伸斷裂應變，γp和γlim為當前時刻剪切塑性應變及剪切斷裂應變，c(t+Δt)和σt(t+Δt)為下一時刻的黏聚力及抗拉強度。

本文采用非關聯流動法則對實體單元的塑性流動行為進行描述，故需在流動法則中引入剪脹角Ψ。

（2）界面本構

邊坡巖體中結構面的密度、產狀（傾向、傾角）及力學參數等均會對邊坡的爆破穩定性產生重要影響，在數值模擬時需要重點考慮。其中，巖體結構面的密度、產狀等幾何參數，在建立邊坡模型及劃分網格時直接進行考慮；對于力學參數，則需要通過結構面力學本構進行體現。

本文采用考慮應變率效應的斷裂能本構，對巖體結構面及巖塊虛擬界面在爆炸載荷下的損傷破裂行為進行精確描述。

首先計算法向試探接觸力，進行拉伸破壞判定，并對法向接觸力及抗拉強度進行修正，為

式中，Fn(t1)為下一時步法向接觸力；Ac為虛擬界面接觸面積；σt0，σt(t0)及σt(t1)為虛擬界面上抗拉強度的初始值、當前值及下一時刻的值；Δun為虛擬界面上法向相對位移的當前值，Gft為虛擬界面拉伸斷裂能（單位：Pa?m）。

接著進行切向試探接觸力的計算，剪切破壞的判定，并對切向接觸力及黏聚力進行修正，為

式中，Fs(t1)為下一時步切向接觸力；?為界面上的內摩擦角；c0，c(t0)及c(t1)為界面上黏聚力的初始值、當前值及下一時刻的值；Δus為界面上切向相對位移的當前值；Gfs為界面的剪切斷裂能（單位：Pa?m）。

當考慮應變率的影響時，初始黏聚力及初始抗拉強度與應變率之間有如下函數關系

式中，γ˙ 及ε˙ 為界面上的剪切應變率及拉伸應變率，e1及e2為界面應變率效應指數，本文取e1=e2=1/3。

1.3 爆生氣體膨脹模型及與圍巖的耦合

本文假定爆轟過程瞬間完成，并借助朗道–斯坦紐科維奇方程（γ 率方程）對爆生氣體的膨脹過程進行刻畫，該方程可表述為

式中，取γ 為3，γ1為4/3，P0和P分別表示初始時刻、當前時刻的爆生氣體壓力，V0和V分別表示初始時刻、當前時刻的爆生氣體體積。Pk和Vk分別為爆生氣體的臨界壓力及臨界體積，Pk的表達式為

式中，Qw為炸藥爆熱（J/kg），ρw為炸藥密度（kg/m3）。

P0的表達式為

式中，D為爆轟速度（m/s）。

爆生氣體與圍巖的耦合采用接觸耦合方式實現，通過法向耦合彈簧，將爆生氣體產生的壓力漸進地傳遞給圍巖；圍巖單元間的法向彈簧及切向彈簧在爆炸壓力的作用下，出現拉伸變形及剪切變形，進而出現拉伸斷裂及剪切斷裂。爆生氣體與圍巖的耦合示意圖如圖2 所示。

圖2 爆生氣體單元與巖石單元的耦合Fig.2 Coupling between blasting gas elements and rock elements

爆生氣體與圍巖的接觸力計算公式可表示為

式中，Fc為法向接觸力，Kn為接觸剛度，Δun為法向接觸位移（拉為正，壓為負）。

由于本文重點關注爆炸應力波對周邊順層邊坡漸進破壞規律的影響，故本節提出的爆生氣體與圍巖的耦合模型僅考慮了爆生氣體對圍巖巖壁的沖擊作用（爆炸應力波產生的主要來源），但未考慮爆生氣體楔入巖體裂縫的過程。

1.4 力學模型的數值化嵌入

采用C++ 語言，將上述邊坡巖體本構模型及爆生氣體膨脹模型代碼化，嵌入至CDEM 的主體程序中。其中，巖體本構模型將拆分成單元本構及接觸面本構，分別從CDEM 程序中的單元本構類及接觸面本構類進行派生；同樣，爆生氣體膨脹模型也將從單元本構類中進行派生，爆生氣體與圍巖之間的接觸模型則從接觸面本構類中進行派生。

力學模型數值化嵌入CDEM 主程序后，將依附于主程序求解框架進行動態求解，單步求解流程如圖3 所示。圖中紅色虛線框部分為本文主要涉及的力學模型及求解流程。

圖3 單步求解流程Fig.3 Solution flow in each time step

2 計算方法及力學模型的案例驗證

對文獻[22]中的爆破漏斗實驗進行數值模擬。該文中的爆破對象為人工材料，由水泥、石膏、河砂、碎石組成，質量分數為1.00∶0.45∶3.82∶3.82。將該人工材料倒入高強度PVC 模具中，在自然條件下28 d 后制成直徑57 cm，高度50 cm的圓柱體試件。在模型中部放入4.5 g TNT 炸藥，埋深17 cm。具體的試驗步驟、方法及試驗材料介紹詳見文獻[22]，此處不再贅述。

為了模擬爆破漏斗試驗，建立了相同尺寸的圓柱體模型（圖4），并用162 744 個四面體單元進行離散。在模型的底部施加全約束邊界條件，在模型的側面施加法向約束條件，模型的頂部為自由面。

圖4 爆破漏斗數值模型Fig.4 Numerical model of blasting crater

人工材料數值模型由單元及虛擬界面組成，單元采用線彈性本構進行描述，虛擬界面采用斷裂能本構進行描述。TNT 炸藥采用朗道爆生氣體膨脹模型進行描述，炸藥參數取值主要來自文獻[22]，其中炸藥密度為1400 kg/m3，爆速為4700 m/s，爆熱為3.0 MJ/kg。炸藥及人工材料間采用法向接觸彈簧進行耦合。人工材料的參數如表1 所示，材料參數主要根據文獻[22]中的材料配合比及文獻[23-24]中關于材料配合比與材料參數的關系綜合確定。

表1 人工材料計算參數Table 1 Computation parameters of artificial material

起爆后，材料的漸進破壞過程及破碎塊體的運動過程如圖5 所示。由圖可以看出，起爆后0.5 ms 時頂面出現了若干徑向裂縫；起爆后3 ms時頂面中部開始隆起；6 ms 后頂部巖體出現大量破裂，內部碎片向外噴出；20 ms 后爆坑已經基本形成。

圖5 人工材料的漸進破壞過程及碎塊運動過程Fig.5 Progressive failure and fragments movement process of artificial material

數值模擬和物理實驗的爆坑對比如圖6 所示。從圖中可以看出，數值模擬得到的爆坑形態、直徑與物理實驗得到的爆坑形狀及直徑基本相同，驗證了CDEM 模擬巖土爆破問題的準確性和合理性。

圖6 數值模擬與物理實驗的爆坑形態對比Fig.6 Comparison of craters between numerical simulation and physical experiment

模型頂部中點的位移時程曲線如圖7 所示，其中實驗部分的位移–時間曲線由高速攝影拍攝的視頻解析所得。由圖可得，數值計算獲得的頂部位移隨時間的變化曲線與現場實驗的結果基本一致，平均誤差在15%左右，滿足工程精度要求。

圖7 模型頂部中點位移時程曲線的對比Fig.7 Comparison of displacement of middle point on the top side

需要說明的是，文獻[22]中的實驗對象為人工均質材料，并沒有設置結構面。本文對該模型實驗進行模擬對比，主要是為了驗證本文提出的巖體損傷破裂模型及爆生氣體–固體耦合算法的正確性。從圖6 和圖7 的對比分析可以看出，本文提出的力學模型及耦合方法用于模擬爆破載荷下巖體的動力學行為是適宜的，計算精度可以保證。

3 爆破開采對順層邊坡穩定性的影響

3.1 計算模型及參數

以鞍千礦業啞巴嶺礦區的臺階邊坡爆破穩定性分析為例。礦區內出露的地層主要為前震旦紀鞍山群及遼河群變質巖系。巖層主要包括赤鐵礦、片巖、灰巖、千枚巖等，地層走向140°，傾向為SW，傾角小于50°。礦區內巖體受結構面切割，形成多處順層臺階邊坡。

概化后的臺階邊坡計算模型及炮孔位置如圖8 所示。邊坡的臺階高度12 m，臺階坡面角65°，平臺寬度2.6 m，邊坡幫坡角55°；炮孔直徑25 cm，孔深15 m，填塞7 m，炮孔到坡腳的距離為6 m。該圖中，編號1，2，3，4，5 表示5 種不同傾角的結構面，結構面傾角分別為10°，20°，30°，40°及50°。

圖8 結構面邊坡數值模型Fig.8 Numerical model of jointed slope

炮孔內炸藥為現場混裝的乳化炸藥，采用爆生氣體膨脹模型進行描述，爆生氣體與圍巖的耦合采用接觸耦合方式進行模擬。乳化炸藥的密度為1150 kg/m3，爆速為4250 m/s，爆熱為3.4 MJ/kg。

計算過程中，為了消除爆炸應力波在人工截斷邊界處的虛假發射現象，本文在模型的底部及四周引入了黏性邊界。

臺階邊坡的力學參數來自文獻[14]，具體如表2 所示。

表2 鐵礦邊坡實體單元計算參數Table 2 Computation parameters of solid elements of iron ore slope

順層結構面的結構面角度為50°，黏聚力為0.1 MPa，抗拉強度為0.1 MPa，內摩擦角為30°，拉伸斷裂能為10.0 Pa?m，剪切斷裂能為50.0 Pa?m。

采用單因素法，分別探討不同爆破距離、不同順層面傾角及不同順層面強度下，爆破對后方臺階邊坡穩定性的影響。討論某一變量的影響時，其他參數取上述基礎參數。

對于結構面傾角及爆破距離，將在建立邊坡幾何模型時直接考慮。其中，結構面建模時，以結構面位置及傾角為界，將邊坡切割為不同的體；炮孔建模時，根據設定距離建立圓柱形炮孔，并與邊坡模型進行求交操作。對于順層面強度，則在數值計算中通過考慮應變率效應的結構面斷裂能本構進行體現。

3.2 爆破距離的影響

研究爆破距離對貫穿性結構面邊坡的影響規律，共探討6 種爆破距離的影響。結構面角度為50°，爆破距離分別為3 m，6 m，9 m，12 m，15 m 和18 m。典型爆破距離下結構面邊坡的水平方向位移云圖如圖9 所示（圖中黑線表示破裂面）。由圖可得，隨著爆破距離的增加，結構面貫穿程度逐漸減小，但變化并不明顯。

圖9 不同爆破距離下結構面邊坡水平方向位移云圖Fig.9 Horizontal displacement contour of jointed slope under different blasting distances

結構面的破裂面積隨爆破距離的變化曲線如圖10 所示。從曲線可看出，當爆破距離從3 m增大至18 m，結構面的破裂面積從1060 m2減小至820 m2，減小比例23%。

3.3 結構面角度的影響

共探討5 種結構面貫穿角度的工況，分別命名為工況1～5（如圖8 所示），工況1～5 對應的貫穿結構面角度分別為10°，20°，30°，40°及50°。

數值計算結果表明，當炮孔距離坡腳為6 m時，爆破載荷對工況1～4 的結構面邊坡無影響，爆破完成后邊坡仍處于穩定狀態，結構面也并未出現破壞。但該爆破距離下的爆破載荷對工況5的結構面邊坡（結構面角度為50°）影響較大，導致該邊坡的結構面出現了整體破壞，并最終造成邊坡失穩。爆破作用下結構面的破壞模式主要為拉–剪復合型破壞，當結構面完全破壞后，上覆巖體將沿著結構面發生整體滑動，如圖11(a)所示（圖中藍色表示未破壞，黃色表示拉伸破壞，紅色表示剪切破壞）。從位移云圖（圖11(b)）可以看出，結構面兩側的位移出現了不連續，表明邊坡已沿著結構面發生了整體滑移，邊坡處于失穩狀態。

圖11 結構面角度50°時的破壞狀態Fig.11 Failure state of structural surface with the dip angle 50 degree

3.4 結構面強度的影響

取炮孔到坡腳的距離為15 m，結構面角度為50°，研究結構面強度對邊坡穩定性的影響。文中采用Mohr–Coulomb 準則描述結構面的壓剪破壞，主要涉及的強度參數包括黏聚力及內摩擦角；采用最大拉應力準則描述結構面的拉伸破壞，主要涉及的強度參數為抗拉強度。本節共選取6 組結構面強度進行探討，每組強度參數中的內摩擦角保持一致，為30°，且每組中結構面的黏聚力和抗拉強度取相同值（即c=σt，以下統稱為結構面強度）。所選取的結構面強度值分別為0.2 MPa，0.6 MPa，1.0 MPa，1.4 MPa，1.8 MPa 和2.2 MPa。

結構面強度為0.2 MPa，0.6 MPa，1.0 MPa和1.4 MPa 時的邊坡水平方向位移云圖如圖12所示。從圖中可以看出，結構面強度對結構面破壞程度的影響較大，隨著結構面強度的增加，結構面貫穿程度逐漸減小。

圖12 不同結構面強度下結構面邊坡水平方向位移云圖Fig.12 Horizontal displacement contour of jointed slope under different structural surface strength

結構面破裂面積隨結構面強度的變化曲線如圖13 所示。從曲線可看出，結構面破裂面積隨著結構面強度的增大而逐漸減??；當結構面強度從0.2 MPa 增大至1.4 MPa，破裂面積從1156 m2減小至337 m2，減小了71%。

圖13 爆破載荷下結構面破裂面積隨結構面強度的變化曲線Fig.13 Relationship between fracture area and strength of structural surface under blasting load

3.5 討論

炮孔內炸藥起爆后，爆炸能量一部分用于臺階自由面附近的巖體破碎，另一部分則向遠處傳播。當爆炸應力波沿著臺階逐漸向后方傳遞，并進一步傳遞至臺階上方邊坡時，若邊坡內存在軟弱結構面，在直達剪切波及坡面反射拉伸波的雙重作用下，或將導致結構面出現拉剪復合型破壞，并進一步誘發順層邊坡出現整體性失穩滑動。

文章借助CDEM 方法，利用單元表征巖體的變形特性，利用單元間的界面表征結構面的斷裂滑移效應，詳細探討了爆破距離、結構面產狀及結構面強度等因素對順層邊坡結構面破壞程度的影響。根據前述計算分析可以看出，由于坡腳處的動應力最大，且坡腳處所受的切向應力最大，因此順層結構面的破裂失穩首先出現在坡腳附近。一旦坡腳處出現破裂，結構面上的整體平衡被打破，裂縫將在動力作用下進一步向坡體中上部的結構面發展。隨著到坡腳距離的增大，動應力峰值逐漸減小，對結構面的擾動也逐漸減小，一旦擾動應力幅值小于結構面上未破裂區域的強度，弱面的破裂失穩過程隨即停止。若減小炮孔到邊坡坡腳的距離，或提升單孔裝藥量，即使爆炸應力波將隨著邊坡高程的增加逐漸減少，但仍然有足夠的能力誘使順層結構面出現整體破裂，屆時順層滑體將會沿著軟弱結構面整體下滑，出現順層滑坡災害。爆破載荷下順層邊坡的失穩機理示意圖如圖14 所示。

圖14 爆破載荷下順層邊坡的失穩機理Fig.14 Failure mechanism of bedding rock slope under blasting load

4 結論

本文以CDEM 為基礎，通過引入考慮應變率效應及應變軟化效應的摩爾–庫倫單元本構模型、考慮應變率效應的界面斷裂能模型、以及爆生氣體絕熱膨脹及與圍巖的耦合模型，實現了炸藥起爆—爆炸應力波傳播—邊坡漸進破裂的全過程模擬。通過爆破漏斗的數值試驗，證明了本文所提出的力學模型及所編制的計算程序在刻畫爆炸載荷下巖體損傷破裂過程方面的優勢。

借助本文提出的力學模型，開展了爆破作用下順層臺階邊坡漸進破壞規律的研究，揭示了爆破載荷下順層邊坡的失穩破壞機理，探討了炮孔到邊坡的距離、結構面產狀、結構面強度等因素對順層邊坡穩定性的影響規律，所得結論如下。

（1）在爆炸應力波的直接壓剪作用及坡面反射波的拉伸作用下，順層結構面易發生拉剪復合型破壞；一旦結構面整體破壞，上覆順層滑體將沿著結構面出現整體性滑動。

（2）結構面傾角越大、結構面強度越低、炮孔到邊坡的距離越小，順層邊坡越容易在爆炸載荷下出現動力失穩；在文章計算參數范圍內，相較于炮孔到邊坡的距離，結構面強度對邊坡爆破穩定性的影響更為顯著。

本文揭示的順層臺階邊坡爆破失穩機理，給出的爆破距離、結構面產狀、結構面強度等參數對順層邊坡爆破穩定性的影響規律，將有助于工程設計人員及現場工作人員確定合理的臺階邊坡參數及爆破安全距離，并將對邊坡爆破穩定性的預警分析起到一定的指導作用。

后續將重點開展反傾邊坡、塊狀巖體邊坡等巖質邊坡在爆炸載荷下的漸進破壞規律分析，通過調整邊坡結構面參數、炸藥參數及孔網參數等，開展大量的仿真分析，形成爆破載荷下邊坡穩定性及變形破壞特征的仿真數據庫，借助機器學習算法構建邊坡爆破分析降階模型，以期為邊坡爆破穩定性的設計評估提供快速分析的工具。