基于CDEM的巖石基坑爆破效果的數值模擬

任慧敏，馮 春，唐昊天，張大帥，王 彪，趙紅華

(1.工業裝備與結構分析國家重點實驗室，遼寧 大連 116024；2.大連理工大學工程力學系，遼寧 大連 116024；3.中國科學院力學研究所流固耦合系統力學重點實驗室，北京 100190；4.中鐵八局集團有限公司城市軌道交通分公司，成都 610000)

隨著基坑開挖深度的不斷增加，開挖面臨的地質條件也更加惡劣，開挖到巖石層的基坑工程大多采用了延時爆破施工[1-4]。延時爆破又稱為毫秒爆破，它在控制地震效應、擴大爆破規模、控制爆破塊度和提高爆破效果方面，以及充分利用爆能、降低藥耗等方面均起著重要作用?？组g毫秒延時間隔的長短是影響延時爆破效果的重要參數，國內外研究學者已經做過很多關于爆破效果相關的研究工作。楊仁樹等[5]采用新型試驗系統，用有機玻璃代替脆性巖石進行了延時爆破裂紋試驗，驗證了延時時間對爆生裂紋擴展有顯著影響。Yan等[6]根據連續介質力學原理和動力非線性有限元理論，推導了炸藥和巖石的控制方程，采用罰函數法結合腐蝕算法，建立了相應的有限元方程，并利用LS - DYNA分析了海底硬巖在延時爆破下的動力響應和爆破動力學特性，結果顯示，合理采用排間延時爆破，可使礦石得到徹底破壞，從而提高巖石爆破質量。Choudhary等[7]利用粒徑為3～7 mm的骨料或篩分鉆屑和合理的行間延時時間，減少鉆鋌產生的巨礫，增加淤泥堆的松散度，從而提高爆破效率及其效果。Wang等[8]基于哈里斯爆破數學模型和疊加法，提出了一種裝藥結構多平面檢測優化算法(MDO)，并通過大直徑深孔爆破作業驗證MDO算法能改善巖石破碎，提高炸藥能量利用率。Yuan[9]等建立1個帶有2個圓形孔的矩形顆粒組合體模型，以相鄰兩孔間距和延時爆破延遲時間為主要變量，分別考慮系數為1.5的不耦合裝藥，發現孤立粒子或約束弱化粒子的碰撞作用對裂隙網絡的形成起著至關重要的作用。潘博等[10]對不同節理角度的試件進行SHPB試驗，并結合核磁共振成像系統得到最利于巖石破碎的荷載與節理的成角范圍。

從以上研究成果可以看出，目前針對巖石基坑爆破效果方面數值模擬研究的成果仍然不多。因此，本文依托大連地鐵 5 號線虎灘公園站基坑爆破工程，利用連續非連續方法CDEM中的塊體動力學軟件(GDEM-BlockDyna)，研究巖石基坑模型在延時間隔時間和炮孔間距2個變量下的爆破效果及規律，并通過應力、振動速度、平均破碎尺寸、極限破碎尺寸、大塊率及系統破裂度等指標進行綜合分析。

1 連續-非連續數值模擬方法及力學模型

1.1 方法概述

GDEM 力學分析軟件，以連續-非連續數值模擬方法(CDEM)為基礎，通過 GPU/CPU的并行加速，可以分析材料在各類靜、動荷載作用下的損傷破裂過程，實現工程尺度下材料的漸進破壞過程。CDEM是一種網格粒子高度融合的連續-非連續顯式數值分析方法。該方法將連續介質數值方法與非連續介質數值方法進行深度融合，在能量層面實現了有限元、離散元及無網格算法的統一。該方法以Lagrange方程為理論基礎，其表達式為

(1)

式中：ui，vi為廣義坐標和廣義速度；L為Lagrange函數，它是系統的能量；Qi為廣義力。

由于是在積分意義下討論連續-非連續計算過程，并沒有考慮離散格式，所以無論是網格還是粒子(見圖1)，只要在其上引入積分點，就可以開展相應的工作。這就是連續-非連續方法囊括了有限元、塊體離散元、顆粒離散元以及無網格等眾多網格粒子這一類統一方法的原因。

圖1 網格及粒子積分點Fig.1 Grid and particle integral points

GDEM塊體動力學軟件(GDEM-BlockDyna)是一款基于多核CPU并行的顯式動力學高效數值仿真軟件，利用JavaScript的各類函數進行求解流程的精確控制。該軟件以CDEM為計算核心，利用塊體表征材料的連續介質特性，利用塊體間的界面表征材料的非連續介質特性，塊體邊界及內部產生斷裂，完成材料漸進破壞過程的仿真過程。該軟件已經在巖土、采礦、爆破、隧道、油氣、水利、地質等多個領域成功應用，服務于相關行業的科研及生產。

1.2 JWL點火爆炸模型

BlockDyna程序描述高能炸藥爆轟產物壓力-體積關系采用JWL狀態方程[11]。高能炸藥爆轟產物的單元壓力p由狀態方程(2)求得，具體可表述為

(2)

式中：V為相對體積；E0為初始內能密度，J/m3；A、B、R1、R2、ω為圓筒試驗確定得出的常數[12-14]。

根據王鴻運[15]的研究成果，本次模擬中JWL狀態方程各參數取值如表1所示。選擇2號巖石乳化炸藥[16]進行爆破模擬，計算所得炸藥參數如表2所示。

表1 炸藥JWL 狀態方程參數

表2 炸藥參數

爆炸模型的計算流程如圖2所示。

圖2 爆炸模型計算流程Fig.2 Calculation procedure of explosion model

1.3 本構模型

采用巖體損傷破裂的力學模型，該模型由單元及虛擬接觸面組成，即單元采用線彈性模型進行描述，輸入包括密度及彈性模量等材料參數；單元間的虛擬接觸面采用斷裂能模型進行描述[17-18]，輸入參數包括黏聚力、內摩擦角、抗拉強度、法向連接剛度、切向連接剛度、拉伸斷裂能、剪切斷裂能。一般情況下，單元及虛擬接觸面上的黏聚力、內摩擦角、抗拉強度的取值一致。

單元線彈性本構模型為

(3)

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

接觸面的損傷斷裂計算首先要確定界面上下一時步的法向及切向試探接觸力：

(4)

式中：Fn、Fs分別為法向、切向接觸力；kn、ks為單位面積上法向、切向接觸剛度，Pa/m；Ac為接觸面的面積；Δdun、Δdus為法向、切向相對位移增量。

采用式(5)判斷是否發生拉伸破壞、修正法向接觸力及抗拉強度，為

(5)

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

采用式(6)判斷是否發生剪切破壞、修正切向接觸力及黏聚力，為

(6)

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

2 數值計算

2.1 工程概況

大連地鐵5號線全長23.8 km，設車站18座，其中虎灘公園站為全線第二站。車站中心里程 YK2+646.963，車站起點里程 K2+577.963，車站終點里程K2+756.958，長178.995 m，寬約 20.0 m，高約18 m，設計底板高程-14.30 m，覆土厚 1.6 ～ 4.8 m，擬采用明挖法施工。

2.2 計算方案

在基坑開挖的過程中，最上層的素填土采用機械開挖方式，素填土的深度為3 m；對于強風化石英巖和中風化石英巖，采用爆破開挖方式，強風化石英巖和中風化石英巖的深度分別為1.5 m和22.5 m，本文著重分析中風化石英巖部分。由于爆破計算所劃分的網格尺寸在厘米量級，若是建立三維數值模型，計算規模及耗時都是巨大的。故為了節約計算時間及成本，建立尺寸為20 m × 10 m的二維平切面模型(見圖3)。該爆破模型內設置3個起爆點，炮孔直徑設為250 mm，藥包間距為L。模型的頂面設置為自由面，炮孔到自由面的距離為5 m，剩余的側面以及底面設置為無反射邊界條件(黏性邊界)，避免爆炸在外邊界產生的應力波進行虛假反射，圖3a中模型上半部紅色邊框范圍內為爆破塊度統計區域(域A)。在BlockDyna中導入ABAQUS軟件劃分的模型網格(見圖3b)，在起爆點周圍進行網格的局部加密。

圖3 爆破模型Fig.3 Blasting model

單純考慮不同延時時間影響時，固定孔距L=5 m，目前國內部分延時爆破工程中多采用Δt=15～75 ms[19-20]，經研究確定了合適的延時爆破時間，間隔時間依次是25、35、50、75 ms。

單純考慮不同炮孔間距影響時[21-22]，固定延時間隔時間為35 ms，炮孔間距依次取為3、4、5、6、7 m。

根據巖土的勘察報告，給出了中風化石英巖物理力學指標參數(見表3)。

表3 巖石參數

2.3 爆破效果評價指標

本文在考慮巖石基坑質點應力時程及峰值振速的同時，為了評定不同爆破參數對巖石的爆破效果，結合爆破后爆破塊度統計區域的塊度大小分布特點進行計算分析，采用4個評價指標[23]。利用腳本文件輸出塊體級配曲線信息，對數據進行整合后得到各指標具體數值[24-27]。

1) 平均破碎尺寸(d50)：塊體級配曲線信息中通過率為50%時相應的特征尺寸；該值與巖石的爆破碎塊尺寸平均值成正比，該值越小表示模型破碎效果越好。

2) 極限破碎尺寸(d90)：塊體級配曲線信息中通過率為90%時相應的特征尺寸；該值與巖石爆破后的大碎塊尺寸成正比，該值越小表示模型破碎效果越好。

3) 系統破裂度(Fr)：數值模型虛擬界面產生破裂的面積占總面積的比值；該值越大表示模型破碎越完整。

4) 大塊率(Br)：塊體級配曲線信息中特征尺寸大于0.9 m的巖石體積之和與模型總體積的比值，該值與模型破碎效果成反比。

3 計算結果分析

巖石基坑模型在基坑穩定線彈性分析的基礎上施加爆炸沖擊荷載，調用BlockDyna開展顯式動力學模擬，設定每100時步向平臺推送計算結果，可得到模型單元位置處位移、速度、加速度、應力等云圖信息。在計算模型2、3炮孔連線中點與頂部自由面的交點處設置監測點，并開展數值計算分析。有關CDEM方法計算爆破的精度及合理性已經被大量文獻驗證[28-29]，此處不再贅述計算的可靠性。

3.1 延時時間的影響

4種延時爆破模型監測點處的最大主應力時程曲線特性趨于一致，選取延時35 ms的應力時程曲線(見圖4)進行表述，可以看出，選取的監測點處所受到的最大正應力正向峰值明顯高于負向峰值，巖體以拉應力破壞為主導。最大主應力峰值分別為7.585 1、7.351 54、5.802 1、5.847 9 MPa，均大于爆破區域巖石的抗拉強度1.4 MPa，能達到巖石破壞的要求。延時時間對爆破前期作用的峰值差別不大，監測點的最大主應力峰值在第2個炮孔起爆附近，爆破后期曲線的躍增逐漸減小。間隔時間為25 ms和35 ms監測點處應力波波形產生重疊，巖石受疊加應力而易于破碎，爆破效果更好。

圖4 延時35 ms的應力時程Fig.4 Stress time-history of 35 ms delay

4種延時爆破模型監測點處y方向(垂向)振動速度時程曲線如圖5所示，在第二個炮孔起爆后振動速度急速增長，測點在y方向速度最大值分別為：17.32、16.67、18.89、17.03 m/s，此后速度基本保持不變。各模型速度最大值分別出現在41、53、67、106 ms，均延遲于第2炮孔起爆時間，說明測點位置處應力波持續沖擊，使巖石塊體加速產生位移。4種延時爆破模型監測點處振動速度最大值在16.67 ～18.89 m/s，相差最大約2.22 m/s，采用35 ms的延時間隔時間質點振速最小。

圖5 y方向不同延時的振動速度Fig.5 Vibration velocity with different delay in y direction

分析監測點的應力曲線、振速曲線可知，延時35 ms時巖石爆破效果較好，因此選取延時35 ms分析巖石爆炸破碎的演化過程(見圖6)。由圖可知，第1個藥包起爆后，沖擊波自爆源向四周傳播，炮孔附近的巖石被擠壓粉碎，裂紋沿著網格呈放射性布滿炮孔周圍；隨著時間的推移，沖擊波衰減為壓應力波，使巖石質點產生徑向位移，在高壓氣體的膨脹、擠壓作用下形成倒錐形裂紋；35 ms時第2個炮孔開始起爆，炮孔間應力波相互疊加，以致大塊巖石被擊碎成細小的粉粒；當時間大于86 ms后，破碎基本形成。

圖6 延時35 ms時不同時刻的總位移云圖Fig.6 Displacement contours at different time with 35 ms delay

延時35 ms爆破時，基坑模型的最終破碎形態如圖7所示。由圖7可知，逐孔起爆激活了幾乎所有的預設虛擬界面，縱橫交錯的裂隙將巖石切割破碎，炮孔上側的巖石破碎效果明顯，下側巖石局部產生拉伸裂隙，爆區頂部地表有大范圍的隆起、位移。由模型爆破的最終形態可以看出，當延時間隔為35 ms時，巖石破碎塊體較致密，爆破效果最好。延時25 ms爆破效果最差，較短的延時間隔不利于應力波在巖石裂隙內的傳播和衰減，導致巖石局部產生過于破碎情況，而其他區域破碎不充分。

圖7 延時35 ms的爆區破碎效果Fig.7 The fragmentation effect in blasting area of 35 ms delay

不同延時爆破的塊體級配曲線如圖8所示。隨著特征尺寸的增大，塊體的通過率逐漸增大；當延時間隔為35、50、75 ms時，通過率達到100%時的特征尺寸依次增大，巖石爆破效果逐漸變差，其中延時25 ms時的爆破效果最差(大塊較多)。將監測數據進行處理得到平均破碎尺寸、極限破碎尺寸、系統破裂度、大塊率等評價指標如表4所示?？梢钥闯?，延時為35 ms時，巖石破碎充分，碎塊尺寸較小，爆破效果最好。

圖8 不同延時間隔的塊體級配Fig.8 Grading of blocks with different delay intervals

表4 不同延時的爆破破碎評價指標取值

不同延時爆破模型爆破塊度統計區域的系統破裂度時程曲線如圖9所示，分析可知：系統破裂度隨著時間的增加而逐漸增大，曲線在每一個炮孔起爆后均有突變且破裂度增加速度加快；爆炸初期系統破裂度增加速度較快，當最后一個炮孔起爆后，破裂度的曲線斜率迅速變緩； 爆破完成后，25、35、50、75 ms延時爆破模型最終的系統破裂度分別為53.7%、56.9%、56.4%、56.0%，延時35 ms的系統破裂度較大，表征爆破效果更佳。

圖9 不同延時下系統破裂度時程Fig.9 Time history of system fracture degree under different delay

3.2 炮孔間距的影響

炮孔間距L=5 m爆破后監測點處的最大主應力時程曲線如圖10所示。不同炮孔間距下，監測點最大主應力分別為5.900 1、8.008 3、8.262 0、6.138 7、 11.957 3 MPa，均超過巖石抗拉強度。在第2個炮孔起爆后的一段時間內，測點處應力波動密集；L=5 m和L=6 m時，監測點處應力波重疊效果明顯，受到的拉伸破壞作用加強，巖石裂紋發育充分。

圖10 炮孔間距5 m的應力時程Fig.10 Stress time-history of 5 m hole spacing

5種炮孔間距監測點處y方向振動速度時程曲線如圖11所示，曲線整體走勢一致呈上升狀態，測點在y方向速度最大值分別為：20.16、19.02、17.82、14.61、11.97 m/s；隨著L的增大，測點處y方向峰值振動速度逐漸減小，相較于L=3 m，相差最大約8.19 m/s。

圖11 y方向不同炮孔間距的振動速度Fig.11 Vibration velocity with different hole spacing in y direction

不同炮孔間距下的基坑模型的最終破碎形態相似，選取L=5 m進行表述(見圖12)。由圖可知，隨著孔距的增加，爆區裂紋擴展更加充分，碎塊尺寸更加均勻。當L=3 m時，由于炮孔間距過小，阻礙了應力波的傳播，導致巖石局部產生過于破碎情況，而其他區域破碎不充分。當L=7 m時，由于炮孔間距過大，爆炸應力破壞效應減弱，導致炮孔周圍出現較密集區域，其他區域碎塊尺寸較大。由模型爆破的最終形態可以得出，當L=5 m和6 m時，巖石破碎塊體較致密，爆破效果最好。

圖12 炮孔間距5 m的爆區破碎效果Fig.12 The fragmentation effect in blasting area of 5 m hole spacing

不同炮孔間距的塊體級配曲線如圖13所示?？梢钥闯鰤K體的通過率隨著特征尺寸的增大逐漸增大；隨著炮孔間距的增大，通過率達到100%時的特征尺寸越小，巖石碎塊越均勻；其中L=6 m時的特征尺寸略小于L=7 m的特征尺寸。將監測數據進行處理得到平均破碎尺寸、極限破碎尺寸、系統破裂度、大塊率等評價指標如表5所示。從表5可看出，當L=6 m時，巖石破碎充分，碎塊尺寸較小且均勻，爆破效果最好。

圖13 不同炮孔間距的塊體級配Fig.13 Grading of blocks with different hole spacing

表5 不同孔距的爆破破碎評價指標取值

不同炮孔間距下系統破裂度時程曲線如圖14所示，分析可知：系統破裂度隨著時間的增加而逐漸增大，不同孔距的突變點大致相同；爆破完成后，5種炮孔間距爆破模型最終的系統破裂度分別為58.7%、59.6%、60.6%、60.7%、60.3%，L=6 m和L=7 m的系統破裂度較大，表征爆破效果更佳。

圖14 不同孔距下系統破裂度時程Fig.14 Time history of system fracture degree under different hole spacing

4 結論

1)對不同延時爆破的計算模型進行爆破應力及振動速度分析，根據模擬結果，測點振動速度最大值在16.67～18.89 m/s，延時25 ms和35 ms均可獲得較好的破碎效果，而延時35 ms引起的峰值振速較小。根據地鐵基坑工程的爆破要求，延時35 ms的模型爆破效果更好。

2)延時35、50、75 ms時，平均破碎尺寸、極限破碎尺寸、大塊率隨著延時間隔的增大逐漸增大，系統破裂度隨著延時的增大逐漸減小。從各指標可以看出，巖石爆破模型的損傷破裂效果隨著孔間延時的增大而逐漸變差，延時25 ms爆破效果最差，較短的延時不利于應力波在巖石裂隙內的傳播和衰減，導致巖石局部產生過于破碎情況，而其他區域破碎不充分。

3)炮孔間距L=5 m和L=6 m時，所選取的監測點所受應力波疊加破壞作用明顯；測點振動速度最大值在11.97～20.16 m/s，測點處y方向振動速度峰值隨著炮孔間距L的增大逐漸減小，故炮孔間距L=6 m的爆破效果優于L=5 m的爆破效果。

4)炮孔間距L=3、4、5、6 m時，平均破碎尺寸、極限破碎尺寸、大塊率隨著炮孔間距的增大逐漸減小，系統破裂度隨著炮孔間距的增大逐漸增大?？梢缘贸?，隨著炮孔間距的增大，破碎塊體的均勻性逐漸提高，但炮孔間距L增加到7 m時，巖石的平均破碎尺寸達到最大，炮孔間距過大使得應力疊加效應減弱。

5)改變炮孔間距相對于改變延時時間大塊率從7.9%到4.7%降低了3.2%，系統破裂度從56.9%到60.7%提高了3.8%，改變炮孔間距的爆破效果優于改變延時時間的爆破效果。