基于連續-非連續單元方法的露天礦三維臺階爆破全過程數值模擬

馮 春，李世海，鄭炳旭，崔曉榮，賈建軍

(1.中國科學院力學研究所流固耦合系統力學重點實驗室，北京 100190；2.廣東宏大爆破股份有限公司，廣東 廣州 510623；3.鞍鋼礦業爆破有限公司，遼寧 鞍山 114046)

爆破開采是露天礦采選總成本控制的首要環節，穿孔爆破的成本較低，僅占整個采選總成本的1/15左右；但是，爆破效果的好壞將直接影響到鏟裝、運輸、破碎等工序的生產效率及能耗。因此，爆破階段通過改變爆破參數，增大礦巖的損傷破碎程度，減少大塊、根底、巖墻等不利因素，將有助于提升后續工序的生產效率。

數值模擬是進行露天礦爆破開采優化設計及爆破效果分析的有效手段。目前，真正可用于爆破破巖全過程模擬的軟件主要包括MBM(mechanistic blasting model)軟件[1-3]、DMC(distinct motion code)軟件[4-8]以及Blo_Up(blast layout optimization using PFC3D)軟件[9-11]。MBM及DMC均來自澳大利亞奧瑞凱(Orica)公司，其中MBM是有限元與塊體離散元相結合的2D爆破效果數值模擬軟件， DMC是基于顆粒離散元的2D/3D爆破效果數值模擬軟件。Blo_Up是HSBM(hybrid stress blast model)項目研究成果的集中體現，由ITASCA公司負責產品研發，并由昆士蘭大學負責軟件模擬結果的驗證。Blo_Up采用非理想爆轟模型描述炸藥的爆轟過程，采用裂隙流動模型描述爆生氣體的楔入破巖過程，采用格子模型描述破碎巖塊的碰撞、飛散及堆積過程。

由于MBM、DMC及Blo_Up等露天礦爆破全過程數值模擬軟件并未公開發售，因此，對于爆破破巖問題的數值模擬，還是主要借助一些商業軟件或一些科研人員自己編制的科研代碼，如ANSYS/LS-DYNA、AUTODYN、UDEC/3DEC、FEM-DEM、DDA、CDEM等。丁希平[12]采用LS-DYNA軟件模擬了炸藥單耗及填塞長度對爆破效果的影響；璩世杰等[13]利用LS-DYNA詳細探討了節理角度對預裂爆破成縫效果的影響規律；Hu等[14]在LS-DYNA中實現了SPH(smoothed particle hydrodynamics)與FEM(finite element method)的耦合計算，開展了預裂爆破及臺階爆破中巖體損傷破裂過程的模擬。謝冰等[15]采用AUTODYN 2D與UDEC(universal distinct element code)相結合的方法，探討了節理組與炮孔連線的夾角對預裂爆破效果的影響；周旺瀟等[16]探討了利用3DEC(three-dimensional discrete element code)模擬爆破破巖過程的可行性，提出了考慮爆破塊度的3DEC塊體人工離散方法。Yan等[17]利用3DEC軟件，通過在炮孔破碎圈外側施加等效的爆破載荷，并通過控制不同區域的網格尺寸，研究了三維條件下爆堆的形成過程。Trivino等[18]通過FEM(finite element method)-DEM(discrete element method)方法及跨孔聲波測試法研究了爆炸應力波及爆生氣體聯合作用下巖體中的裂紋擴展過程。甯尤軍等[19]通過在DDA(discontinuous deformation analysis)中引入爆生氣體膨脹模型，實現了巖石爆破過程中炮孔擴張、巖體破壞、塊體拋擲和爆堆形成過程的模擬。鄭炳旭等[20]利用CDEM(continuum-based discrete element method)軟件探討了炸藥單耗對爆破塊度分布曲線及系統破裂度的影響規律。

總體而言，學者們對二維情況下巖石的爆破破碎、飛散及堆積過程已進行了深入研究，但三維情況下的研究依然較少。為此，本文中以CDEM方法為基礎，通過引入爆破破巖相關的算法及模型，開展工程尺度下露天礦三維爆破全過程的數值模擬。

1 數值方法及力學模型

1.1 連續-非連續單元方法概述

連續-非連續單元方法(continuum discontinuum element method，CDEM)[21-22]是一種有限元與離散元相互耦合的顯式動力學數值分析方法。CDEM的理論基礎是拉格朗日方程，為：

(1)

式中：ui、vi為廣義坐標，L為拉格朗日系統的能量，Qi為非保守力所做的功。

為了表征多裂紋的萌生、擴展及交匯貫通過程，在CDEM中引入了虛擬裂縫的概念。虛擬裂縫位于每個有限元單元的邊界上，在斷裂發生之前，通過引入法向罰彈簧及切向罰彈簧可連接兩側的實體單元，并進行力學信息的傳遞；通過在罰彈簧上設置斷裂準則及對應的強度參數，可在虛擬界面上實現拉伸斷裂過程及剪切斷裂過程；斷裂發生后，虛擬界面即轉化為真實的接觸界面，通過賦予相應的接觸模型及接觸參數，即可對接觸面的力學行為進行準確刻畫。

CDEM中的核心控制方程為：

(2)

CDEM采用基于增量方式的顯式歐拉前差法進行問題的求解，整個計算過程中通過不平衡率表征系統受力的平衡程度。求解控制方程(2)，共分為3個步驟：第1步循環每個有限元單元，計算單元的變形力及阻尼力；第2步循環每個接觸面，計算接觸面上的連接力及阻尼力；第3步循環所有節點，計算每個節點的合外力、加速度、速度及位移。

1.2 爆源模型

本文的爆源模型主要采用朗道爆炸模型，該模型的輸入參數包括裝藥密度，炸藥爆速、爆熱及點火點位置。該模型主要基于朗道-斯坦紐科維奇公式(γ率方程)，為：

(3)

式中：γ= 3，γ1= 4/3，p、V分別為高壓氣球的瞬態壓力和體積，p0、V0分別為高壓氣球初始時刻的壓力和藥包的體積，pk、Vk分別為高壓氣球在兩段絕熱過程邊界上的壓力和體積。pk的表達式為：

(4)

式中：Qw為單位質量炸藥爆熱(J/kg)，ρw為裝藥密度(kg/m3)。p0的表達式為：

(5)

式中：D為爆轟速度(m/s)。

本文采用到時起爆的方式模擬點火過程及爆轟波在炸藥內的傳播過程。設某一炸藥單元的重心到點火點的距離為d，炸藥的爆速為D，則點火時間t1=d/D。當爆炸時間t>t1時，根據式(3)對該單元進行爆炸壓力的計算。

鉆孔爆破過程中，隨著炮孔的起爆，爆生氣體逐漸膨脹，并推動巖體做功。巖體在爆炸壓力的作用下逐漸出現裂縫的萌生、擴展，并最終出現貫通自由面的裂縫。一旦出現貫通性的裂縫，炮孔內的氣體將會從裂縫中快速溢出，并導致炮孔內的壓力急劇下降。由于直接模擬爆生氣體在巖體內的流動及溢出過程較為復雜，本文中采用設置炸藥作用時間的方式進行等效模擬。當某炸藥單元起爆后經歷的時間大于炸藥作用時間時，該炸藥單元隨即失效，失效后炸藥單元中的氣體壓力為零，且不再進行爆炸壓力的計算。

1.3 巖體本構

本文中采用的巖體本構為彈性-損傷-斷裂本構，其中在每個有限元單元上施加線彈性本構，輸入的參數包括密度及彈性參數(彈性模量、泊松比)；在虛擬界面上施加損傷-斷裂本構，輸入的參數包括法向連接剛度、切向連接剛度、黏聚力、內摩擦角、抗拉強度、剪切斷裂能及拉伸斷裂能。

利用增量法表述的單元線彈性本構為：

(6)

式中：σij為應力張量，Δσij為增量應力張量，Δεij為增量應變張量，Δθ為體應變增量，K為體積模量，G為剪切模量，δij為Kronecker記號，t1為下一時步，t0為當前時步。

虛擬界面上采用考慮斷裂能的拉剪復合本構進行損傷斷裂的計算。首先采用增量法計算虛擬界面處下一時步的法向及切向試探接觸力，為：

(7)

式中：Fn、Fs為罰彈簧上法向及切向的連接力，kn、ks為單位面積上法向、切向連接剛度(Pa/m)，Ac為虛擬界面的面積，Δdun、Δdus為法向、切向相對位移增量。

采用式(8)進行拉伸斷裂的判斷、法向連接力及拉伸強度的修正,如果:

-Fn(t1)≥σt(t0)Ac

那么：

(8)

式中：σt0、σt(t0)及σt(t1)為初始時刻、本時刻及下一時刻虛擬界面上的拉伸強度，Δun為當前時刻虛擬界面上的法向相對位移，Gft為拉伸斷裂能(Pa·m)。

采用式(9)進行剪切斷裂的判斷、切向連接力及黏聚力的修正，如果:

Fs(t1)≥Fn(t1)tanφ+c(t0)Ac

那么：

(9)

式中:φ為虛擬界面的內摩擦角，c0、c(t0)及c(t1)為初始時刻、本時刻及下一時刻虛擬界面上的黏聚力，Δus為當前時刻虛擬界面上的切向相對位移，Gfs為剪切斷裂能(Pa·m)。

基于式(8)及式(9)，可繪制出虛擬界面上法向及切向的本構曲線，具體如圖1所示。

圖1 虛擬界面上的本構曲線Fig.1 Constitutive curves at virtual interface

1.4 接觸碰撞算法

爆破后，破碎巖塊將以一定的初速度向臨空面拋擲，進而出現飛散、碰撞、堆積等爆破現象。本文中采用半彈簧-半棱聯合接觸碰撞模型[23-24]對大量破碎巖塊的接觸碰撞過程進行快速、精確地模擬。

半彈簧由單元頂點縮進至各棱(二維)或各面(三維)內形成；半棱僅存在于三維情況，由各面面內相鄰的兩根半彈簧連接而成(見圖2)。建立半彈簧時，縮進距離一般取頂點到各棱或各面中心距離的1%～5%(本文中取5%)。由于只有在半彈簧找到目標面、半棱找到目標棱以后，方能形成完整的接觸對并計算接觸力，因此稱之為“半”彈簧及“半”棱(見圖3)。

圖2 半彈簧-半棱示意圖Fig.2 Semi-spring and semi-edge schematics

圖3 兩類接觸對Fig.3 Two types of contact pairs

由于采用了縮進策略，使得半彈簧及半棱均位于面內(二維時位于棱內)，因此半彈簧及半棱均具有各自的特征面積(二維時取單位厚度)，為:

As=Af/Nv,Ae=As,i+As,j

(10)

式中:As及Ae分別為半彈簧及半棱的特征面積，Af為半彈簧及半棱所在母面的面積，Nv為所在母面的頂點個數，As,i及As,j分別為組成半棱的2根半彈簧的面積。

半彈簧-半棱聯合接觸碰撞模型將二維情況下點-點、點-線、線-線3類接觸關系統一為半彈簧-目標棱這一類關系，將三維情況下的點-點、點-線、點-面、線-線、線-面、面-面6類接觸關系統一為半彈簧-目標面及半棱-目標棱這2類關系，從而簡化了計算，提升了接觸檢索效率。

一旦接觸對創建完畢，即可在接觸對上引入傳統的接觸本構，執行接觸力的計算。本文中，當2個塊體處于壓縮接觸狀態時執行線彈性計算，當2個塊體處于剪切狀態時執行庫倫摩擦計算；一旦2個塊體處于拉伸狀態，達到其強度極限時立即斷開接觸彈簧。

2 單自由面爆破實驗的數值驗證

以大理巖爆破模型實驗[25]作為數值模擬對象。該實驗中大理巖塊被加工成尺寸為25 cm×25 cm×25 cm的立方體，并在試樣一側鉆出直徑為0.6 cm、深度為10 cm、間距為5 cm的2個平行鉆孔，2個鉆孔到巖樣頂部自由面的距離均為5 cm(見圖4)。采用DDNP炸藥進行爆破實驗，每孔裝藥0.3 g。由于試件的尺寸較小，為了減少和消除試件四周自由面對爆破效果的影響，在巖樣四周和底部涂上一層黃油，然后用5塊1 cm厚的鐵板夾制(布孔一側不用鐵板夾制)，并用螺栓固定。

圖4 試樣及炮孔尺寸Fig.4 The size of specimen and bore hole

數值計算中采用的巖石密度、彈性模量、泊松比、抗拉強度與文獻[25]中的一致，分別為2 730 kg/m3、61.4 GPa、0.27、4.72 MPa；巖石黏聚力由文獻[25]中給出的單軸抗壓強度(90 MPa)及內摩擦角經驗值(40°)通過Mohr-Coulomb公式計算獲得，為21 MPa；巖石的內摩擦角，拉伸斷裂能，剪切斷裂能，單位面積法、切向剛度，斷裂后滑動摩擦因數等參數文獻[25]中未提供，本文中根據相關資料選取了經驗值，分別為40°、50 Pa·m、100 Pa·m、2×1014Pa·m、0.25。數值計算中采用的炸藥密度、爆速與文獻[25]中的一致，分別為1 000 kg/m3、4 950 m/s；炸藥的爆熱、炸藥作用時間等參數文獻[25]中未提供，本文中根據相關資料選取了經驗值，分別為1.4 MJ/kg、5 ms。

文獻[25]中根據耦合及填塞的不同，共開展了3組實驗，分別為耦合無填塞、耦合填塞及不耦合填塞(不耦合因數為1.6)，本文中重點對比耦合填塞的實驗結果。

建立棱長為25 cm的立方體數值模型，在模型的相應位置設置炮孔，采用15.5萬個四面體單元進行剖分。為了模擬實驗過程中的鐵板夾制作用，將模型底部及四周側面設置為全約束，模型頂部表面自由。

雙孔同時爆破后，不同時刻巖體的破碎情況如圖5所示。由圖5可得，自由表面處的巖體在爆炸應力波及爆生氣體的聯合作用下，破裂為若干大塊；在內部氣體的推動作用下，大塊逐漸向外拱出，并發生解體破碎，而爆源附近的大量碎塊隨之涌出。

圖5 不同時刻巖體的破碎運動情況Fig.5 Fracture and movement of rock at different times

爆破7.50 ms后，對脫離母體的碎塊進行刪除，可獲得爆破漏斗的空間形態，如圖6所示。對爆破漏斗的形態進行測量，可得A-A′剖面中，爆破漏斗的深度為5.8 cm，爆破漏斗直徑為23 cm；B-B′剖面中，爆破漏斗的深度為6.9 cm，爆破漏斗的直徑為16 cm。

圖6 爆破漏斗的形態Fig.5 Shape of crater

圖7 塊度分布曲線Fig.7 Block distributing curves

對爆破后各碎塊的體積及特征尺寸進行統計，繪制出塊度分布曲線對比圖(見圖7)。其中，文獻[25]中的特征尺寸為篩孔直徑，數值計算中各碎塊的特征尺寸Lc的計算公式為：

(11)

式中：Vb為某碎塊的體積，Lmax為該碎塊中各頂點間的最大距離。

由圖7可得，在對數坐標系下，2條曲線的變化規律基本一致，數值計算給出的塊度分布曲線基本呈反“S”型，而實驗給出的曲線呈拋物線型。當通過率為40%～50%時，2條曲線對應的特征尺寸基本一致，為33～56 mm；但數值計算獲得的特征尺寸在20 mm以下及70 mm以上的碎塊比例明顯高于實驗值。

數值計算及模型實驗獲得的碎塊總體積、K50、K80等指標的對比如表1所示。其中，K50、K80分別為通過率為50%及80%時對應的碎塊特征尺寸。由表1可得，數值計算獲得的碎塊總體積(爆破漏斗體積)及K50與實驗值基本一致，誤差為16%以下；K80的數值解與實驗值差別較大，誤差約為52%。

表1 關鍵指標對比Table 1 Comparison of key indexes

3 露天鐵礦臺階爆破模擬

以鞍千露天鐵礦的臺階邊坡幾何尺寸、巖體性質及孔網參數為基礎，建立如圖8所示的3排7列(共21個炮孔)三自由面臺階爆破模型。臺階高度為12 m，臺階坡角為90°(直立邊坡)，炮孔的間排距及抵抗線均為6.5 m，炮孔直徑為25 cm，深度為15 m，堵塞長度為6.5 m。采用逐孔起爆方式，孔底起爆，孔間延時為42 ms，排間延時為65 ms。

圖8 含21炮孔的三自由面臺階爆破模型Fig.8 The bench blasting model with three free surfaces and twenty-one bore holes

對鐵礦石塊體采用線彈性模型進行描述，其密度為3 200 kg/m3，彈性模量為60 GPa，泊松比為0.25。綜合考慮了塊體間界面爆區內既有裂隙的強度，采用考慮強度軟化效應的損傷斷裂模型進行描述，單位面積上的法向及切向剛度為200 GPa/m，黏聚力為12 MPa，內摩擦角為30°，抗拉強度為4 MPa，剪切斷裂能為150 Pa·m，拉伸斷裂能為50 Pa·m，斷裂后滑動摩擦因數為0.58。

爆破所用的炸藥為現場混裝的乳化炸藥，對其采用朗道爆炸模型進行描述，裝藥密度為1 150 kg/m3，爆速為5 600 m/s，爆熱為3.4 MJ/kg，炸藥作用時間為35 ms。

數值計算共分為2個階段：第1階段為靜力平衡階段，采用虛擬質量法獲得模型在重力作用下的靜態應力場，在該階段，模型的底部及四周為法向約束邊界，重力方向豎直向下，局部阻尼系數取0.8；第2階段為爆破破巖階段，在模型底部及四周施加無反射邊界，計算時步取10 μs，局部阻尼因數取0.03，按照預設的起爆順序及延時進行點火起爆，獲得鐵礦石的破碎、拋擲過程及最終的堆積過程。

采用四面體網格對上述模型進行剖分，共剖分網格51.1萬。受計算量限制，底部平臺未進行延伸，寬度僅為8.0～8.5 m，因此爆破后碎塊的運動范圍將超出底部平臺；為了對超出底部平臺的塊體提供支撐，特設置與底部平臺相同高度的剛性面。

爆堆的演化過程如圖9所示。由圖9可得，2個側壁臨空面相交的區域首先發生了臌脹及拋擲，拋擲方向與等時線方向基本一致。隨后，在2個側壁臨空面上分別出現了“爆花”，臺階邊坡的中部明顯鼓出。當時間大于4.370 s后，爆堆基本形成。

沿著第2列、第4列及第6列炮孔的中心對爆堆進行剖切，觀察不同位置處爆堆的表面及內部形態，如圖10所示。整個爆堆基本呈上陡下緩的斜坡狀，底部坡角為15°～30°，頂部坡角為40°～50°；原爆區頂部出現了較大范圍的隆起，最大高度可達2.6 m。第1排炮孔前側的破碎塊體堆積得較松散，后側區域的破碎塊體堆積得較致密；爆區內的巖體在爆炸載荷作用下破碎明顯，被縱橫交錯的裂縫切割為大量碎塊；爆區外側的巖體破碎較輕微，僅在局部區域出現了沿著徑向的拉伸裂縫。

圖9 不同時刻的總位移云圖Fig.9 Displacement magnitude contours at different times

圖10 爆堆剖視圖Fig.10 Section views of muckpile

圖11 破裂度時程曲線Fig.11 History of fracture degree

定義破裂度為數值計算中發生破裂的虛擬界面面積與總虛擬界面面積的比，則破裂度隨時間的變化如圖11所示。由圖11可得，隨著爆炸時間的增長，破裂度迅速增大；當爆炸時間大于0.39 s后(最后一個炮孔的起爆時間為0.382 s)，破裂度的增大趨勢迅速變緩；爆破完成后，模型的最終破裂度約為85.4%。

采用UniStrong手持式GPS(亞米級精度)，在鞍千礦開展大量的爆堆形態測試。圖12給出了南采區典型臺階爆破后的爆堆形態。由圖12可得，除后緣拉裂槽外，數值計算給出的爆后斜坡形態、頂部鼓起等現象與現場的測試結果基本一致，證明利用CDEM開展三維露天臺階爆破全過程模擬可行。

圖12 南采區的典型爆堆Fig.12 Typical muckpiles in south region

4 結 論

(1)在連續-非連續單元方法(continuum-discontinuum element method,CDEM)中引入了朗道爆炸模型、巖體彈性-損傷-破裂模型及半彈簧-半棱接觸碰撞模型，實現了爆炸作用力的精確計算，巖體損傷破裂過程的準確描述以及爆破碎塊群飛散、碰撞、堆積過程的快速分析。

(2)開展了單自由面爆破實驗的數值對比分析，給出了雙炮孔同時爆破作用下，自由表面處巖體的破裂飛散過程，數值計算給出的爆破漏斗體積及塊度分布曲線與文獻中的實驗結果基本一致。

(3)以鞍千礦南采區典型的鐵礦臺階爆破開采為背景，模擬了3排21個炮孔逐孔起爆后，爆區內巖體的損傷破裂過程、破碎塊體的飛散過程及爆堆的形成過程，數值計算給出的爆堆形態與現場的測試結果具有一定的相似性。