黃永輝，孫博，張智宇,3，錢正乾，王軍，李洪超

（1.昆明理工大學(xué) 電力工程學(xué)院，云南，昆明 650500；2.昆明理工大學(xué) 國土資源工程學(xué)院，云南，昆明 650093；3.昆明理工大學(xué) 云南省中-德藍色礦山與特殊地下空間開發(fā)利用重點實驗室，云南，昆明 650093；4.玉溪大紅山礦業(yè)有限公司，云南，玉溪 650302；5.昆明理工大學(xué) 城市學(xué)院，云南，昆明 650051）

隨著科技發(fā)展，爆破理論及技術(shù)正在向智能化、精細化方向發(fā)展，數(shù)值仿真技術(shù)作為重要研究手段之一，其關(guān)鍵點在于算法和材料本構(gòu)，適用于爆破破巖數(shù)值計算的算法主要有有限元（FEM）、流固耦合（ALE）、粒子流（SPH）等，材料本構(gòu)方面主要有TCK、HJC、Yang 及RHT 等.算法和材料本構(gòu)方程不僅決定計算效率還影響結(jié)果的準度，相關(guān)研究一直是學(xué)術(shù)界公認的研究熱點.

近年來大量學(xué)者[1?8]對巖石損傷本構(gòu)模型展開了較為系統(tǒng)研究；張若棋等[9]利用AUTODYN 數(shù)值軟件中混凝土HJC、RHT 本構(gòu)進行失效強度參數(shù)分析計算，提出利用混凝土的特征強度確定失效強度參數(shù)的方法并進行校驗；王秀麗等[10]采用SPH 法與RHT 本構(gòu)得到楔擋分流結(jié)構(gòu)在內(nèi)部爆炸作用下的爆坑形態(tài)特征、損傷區(qū)范圍等；張德生等[11]采用基于光滑質(zhì)點流體動力學(xué)與RHT 混凝土本構(gòu)的方法得出破碎大塊煤的有效方法；PRAKASH 等[12]采用改進后的RHT 本構(gòu)模型對鋼纖維增強水泥基復(fù)合材料板（SFRCC）在沖擊過程中的各種特性現(xiàn)象得出最佳纖維體積與厚度下的沖擊能力峰值；WANG 等[13]提出一種確定RHT 材料模型參數(shù)的方法；李洪超等[14]通過正交試驗與敏感度分析對RHT 本構(gòu)中參數(shù)進行分析.許多學(xué)者[15?22]對爆破荷載下巖石損傷分區(qū)判定進行了研究，盧文波等[23?25]利用RHT 損傷本構(gòu)進行數(shù)值模擬，并輔以相應(yīng)的現(xiàn)場爆破試驗，對深部隧道開挖爆破下的巖石損傷孕育機理進行較全面研究；潘城等[26]通過SHPB實驗與室內(nèi)爆破試驗對光面爆破參數(shù)進行優(yōu)化，有效控制了超欠挖現(xiàn)象.宋肖龍等[27]利用地質(zhì)雷達信號具有短時非平穩(wěn)的特點，通過HHT 法對信號進行去噪處理以提取有效反映損傷特征的瞬時參量，從而獲得圍巖損傷圖像來進行損傷分區(qū).劉閩龍等[28]建立各向異性動態(tài)損傷本構(gòu)進行隧道爆破損傷影響數(shù)值模擬，并基于聲波測試原理對隧道圍巖損傷進行測量，以驗證各向異性動態(tài)損傷本構(gòu)的準確性.賈海鵬等[29]建立了損傷敏感區(qū)間計算模型，結(jié)合實際工況計算出了敏感區(qū)間的數(shù)值，并校核了區(qū)間內(nèi)外的安全藥量.

史卜濤等[30]基于廣義插值物質(zhì)點法提出一種物質(zhì)點強度折減法，為邊坡的穩(wěn)定性分析提供了新的分析思路；張忠等[31]將點火增長方程與物質(zhì)點法相結(jié)合，對多種材料在撞擊屏蔽炸藥方面進行模擬，驗證了物質(zhì)點法在沖擊起爆問題的可行性；張芮瑜等[32]針對強夯作用下的土體變形問題，將應(yīng)力密度相關(guān)土體本構(gòu)模型與物質(zhì)點法相結(jié)合，分析總結(jié)出強夯過程中的能量轉(zhuǎn)化規(guī)律，為相關(guān)問題的研究提供新視角；王宇新等[33]基于沖擊動力學(xué)與爆炸焊接理論，利用物質(zhì)點法對界面波進行數(shù)值模擬并進一步分析研究了界面波形成機理.張智宇等[34]利用高速攝影儀對爆破破巖的物質(zhì)點進行軌跡抓拍，得到了不同時刻的鼓包運動輪廓.

物質(zhì)點法在大變形、高速碰撞等方面的模擬效果表現(xiàn)良好，但目前對現(xiàn)場爆破效果分析方面的應(yīng)用相對較少；RHT 損傷本構(gòu)在爆炸荷載下對于材料的損傷狀態(tài)描述方面較為突出，然而巖石損傷的判定準則主要還是基于《水工建筑體地下開挖工程施工技術(shù)規(guī)范》[35]中的判定標準，判定方法需要進一步優(yōu)化.本文將從RHT 本構(gòu)的損傷定義出發(fā)，結(jié)合巖石損傷破裂理論與損傷判定準則，對變鈉質(zhì)熔巖在爆炸荷載下的巖石碎裂范圍進行推導(dǎo)，選取物質(zhì)點法（material point method）的數(shù)值算法，通過數(shù)值模擬與現(xiàn)場爆破漏斗試驗相結(jié)合的方式驗證了判定范圍的準確性.

1 爆破碎裂區(qū)判定方法優(yōu)化及應(yīng)用

1.1 判定方法優(yōu)化

巖石在爆破荷載作用過程中，炮孔壁附近的巖石受炸藥產(chǎn)生的爆轟波影響，一定范圍內(nèi)的巖石被粉碎，此時粉碎區(qū)的巖石損傷變量Dcrush=1，而粉碎區(qū)之外的巖石，由于爆炸能量的衰減，其造成的損傷逐漸減小，因此根據(jù)損傷程度分為巖石碎裂區(qū)、損傷擾動區(qū)及原巖振動區(qū)，不同破壞分區(qū)示意圖如圖1 所示.

圖1 破壞分區(qū)示意圖Fig.1 Failure zone diagram

目前對于巖石碎裂區(qū)與巖石損傷區(qū)之間的損傷臨界值相關(guān)方面研究存在不足，根據(jù)RHT 內(nèi)嵌公式中可發(fā)現(xiàn)，本構(gòu)模型中關(guān)于損傷參量的描述僅有材料初始損傷值D1及完全損傷值D2，并沒有對于巖石臨界損傷與碎裂區(qū)臨界損傷的判定標準.為建立RHT 本構(gòu)中巖石臨界損傷參量與本構(gòu)模型中參數(shù)之間的數(shù)學(xué)關(guān)系，參考吳政等[36]根據(jù)現(xiàn)有對于巖石損傷變量D的研究成果結(jié)合巖石材料損傷失效準則，對巖石材料損傷失效行為進行研究，將巖石材料其峰值應(yīng)力強度所對應(yīng)的損傷值視為巖石臨界損傷參量Dcr，而材料破壞主要由材料發(fā)生塑性變形導(dǎo)致[37]，因此可由式（1）來表示巖石塑性應(yīng)變與巖石臨界損傷參量的關(guān)系：

RHT 本構(gòu)模型中對于損傷參量D定義為

其中：

對式（1）及RHT 本構(gòu)函數(shù)進行整理，分別選取巖石處于開始壓碎時與完全壓實時的應(yīng)力狀態(tài)分別作為巖石損傷區(qū)與巖石碎裂區(qū)臨界閾值評判標準，可得式（1）中的塑性應(yīng)變與極限應(yīng)變?yōu)?/p>

式（1）～（6）中：εp為 塑性應(yīng)變；εmax為 峰值應(yīng)變；為失效時的塑性應(yīng)變；?εp為失效塑性應(yīng)變和當(dāng)前塑性

1.2 變鈉質(zhì)熔巖力學(xué)試驗及RHT 本構(gòu)參數(shù)

研究以大紅山鐵礦地下爆破工程為背景，試驗區(qū)域巖石為變鈉質(zhì)熔巖，取樣開展室內(nèi)力學(xué)性能測試部分試驗過程見圖2，獲得了主要的物理力學(xué)參數(shù)，根據(jù)RHT 本構(gòu)參數(shù)確定方法的相關(guān)文獻[38]，獲得了變鈉質(zhì)熔巖參數(shù)詳見表1.

表1 變鈉質(zhì)熔巖RHT 本構(gòu)參數(shù)取值匯總表Tab.1 Summarization table of RHT constitutive parameters of sodium metamorphosed lava

圖2 室內(nèi)巖石力學(xué)試驗過程Fig.2 Laboratory rock mechanics test process

1.3 變鈉質(zhì)熔巖碎裂區(qū)判定

根據(jù)學(xué)者的研究成果，巖石損傷范圍與巖石初始損傷、單孔裝藥量、爆心距及初始傳播頻率之間呈復(fù)雜的函數(shù)關(guān)系；參考唐紅梅等[39]總結(jié)得出的地下工程施工爆破過程中圍巖損傷分區(qū)中爆心距r與損傷變量D間的關(guān)系曲線.

式中：a2和f0分別代表衰減系數(shù)、爆破初始頻率.

由RHT 本構(gòu)參數(shù)和式（5）、（6）可求得巖石損傷臨界閾值Dcr=0.11，巖石碎裂臨界閾值Dcf=0.51；由式（7）衰減系數(shù)、爆破初始頻率取值參考熊海華等[40]研究成果分別取9.5×10?4、47 Hz，獲得了RHT 本構(gòu)爆破損傷判定范圍如表2 所示.

表2 RHT 本構(gòu)爆破損傷判定范圍Tab.2 RHT constitutive blasting damage determination range

2 基于物質(zhì)點法爆破碎裂區(qū)數(shù)值模擬驗證

2.1 物質(zhì)點算法理論

物質(zhì)點法[41]通過流體隱式質(zhì)點法將本構(gòu)方程改為在各質(zhì)點上計算的方式，并結(jié)合等效積分弱形式，采用質(zhì)點離散建立動量方程的離散格式，因此稱為物質(zhì)點法.該方法[42]仍采用拉格朗日質(zhì)點和歐拉網(wǎng)格進行雙重描述，它將連續(xù)體離散成一組質(zhì)點，每個質(zhì)點代表一塊材料區(qū)域并擁有該區(qū)域的質(zhì)量、速度、應(yīng)力和應(yīng)變等物質(zhì)信息，作為一種完全的拉格朗日質(zhì)點類方法，在每步中質(zhì)點和計算網(wǎng)格沒有進行相對運動，避免了歐拉法在非線性對流向產(chǎn)生的數(shù)值困難，且物質(zhì)界面易跟蹤，因此物質(zhì)點法在沖擊、流固耦合等設(shè)計大變形或材料破壞方面的數(shù)值計算中效率明顯優(yōu)于有限元法與SPH 法，物質(zhì)點法示意圖如圖3 所示.

圖3 物質(zhì)點法示意圖Fig.3 Material point method diagram

物質(zhì)點法無需考慮熱交換，僅通過求解動量方程即可確定材料狀態(tài)，物質(zhì)點法中連續(xù)體的動量方程為

同時其邊界控制條件與密度離散化方程為

背景網(wǎng)格節(jié)點的運動方程為

式（8）～（11）中：ρ為材料密度；u為位移；σij為Cauchy應(yīng)力；bi為單位質(zhì)量下的體積力；δ為Dirac delta 函數(shù)；mp、xip分別為物質(zhì)點p的質(zhì)量與坐標.

而諸多學(xué)者[43?44]對物質(zhì)點法進行了大量研究與擴展，證明了物質(zhì)點法在處理大變形問題方面上具有極大的優(yōu)勢，本文使用搭載物質(zhì)點法和RHT 損傷本構(gòu)的Peneblast 軟件作為后續(xù)爆破漏斗數(shù)值模擬分析算法.

2.2 模擬方案

根據(jù)經(jīng)驗以及現(xiàn)場試驗條件開展5 組變埋深爆破漏斗數(shù)值模擬，設(shè)計參數(shù)見表3.

表3 變埋深爆破漏斗試驗參數(shù)Tab.3 Parameters of variable depth blasting funnel test

2.3 模型及參數(shù)

數(shù)值仿真軟件采用由Fortran90 等語言研發(fā)的三維顯式物質(zhì)點法數(shù)值仿真軟件PeneBlast[45]，該軟件能較充分體現(xiàn)物質(zhì)點法在處理大變形方面模擬的優(yōu)異性，同時支持RHT 本構(gòu)和炸藥本構(gòu).構(gòu)建等比例1/4 對稱模型，具體尺寸為1.6 m×1.6 m×2.0 m，對稱面采用對稱約束，除自由面外的其他邊界采用無反射約束.模型見圖4.巖石本構(gòu)參數(shù)見表1，炸藥采用2 號巖石乳化炸藥，其本構(gòu)參數(shù)見表4.

表4 2#巖石乳化炸藥參數(shù)Tab.4 Parameters of 2# rock emulsion explosive

圖4 數(shù)值模型Fig.4 Numerical model

2.4 模擬結(jié)果分析

（1）爆破漏斗成型過程.

數(shù)值計算完成后對等效塑性應(yīng)變進行輸出，爆破漏斗動態(tài)成型過程如圖5 所示.由圖5 可以看出，在t=3.75～60 ms 時，粉碎區(qū)初步形成，粉碎區(qū)的平均等效塑性應(yīng)變約為2，以壓縮破壞為主；當(dāng)t=90～120 ms時，應(yīng)力波到達自由面并在自由面處發(fā)生反射拉伸波，粉碎區(qū)基本成型，巖石碎裂區(qū)與自由面被拋出部分的巖石等效塑性應(yīng)變高達6.4；t=150 ms 時，爆破漏斗中各損傷分區(qū)基本成型，根據(jù)等效塑性應(yīng)變的云圖顯示，碎裂區(qū)巖石的等效塑性應(yīng)變平均為6.2.

圖5 爆破漏斗動態(tài)成型過程Fig.5 Dynamic forming process of blasting funnel

（2）爆破漏斗破碎分區(qū)規(guī)律.

將被拋擲的物質(zhì)點隱藏處理，并根據(jù)文中推導(dǎo)出的判定范圍對碎裂區(qū)與損傷區(qū)進行范圍標定，詳見表5 中，各組爆破漏斗數(shù)值結(jié)果與損傷范圍分析曲線如圖6 所示.

表5 爆破漏斗數(shù)值模擬數(shù)據(jù)Tab.5 Numerical simulation data of blasting funnel

圖6 爆破漏斗數(shù)值結(jié)果Fig.6 Numerical results of blasting funnel

由表5 和圖6 可知：粉碎區(qū)、碎裂區(qū)、損傷區(qū)平均半徑分別為502 mm、830 mm、1 182 mm，均在理論半徑值所劃分范圍內(nèi).隨著埋深的增加，各分區(qū)半徑呈現(xiàn)出先增大后減小的變化趨勢，分區(qū)半徑在埋深為590 mm 時達到最大，而后半徑中等幅度減??；粉碎區(qū)深度相比炮孔深度多出7.5%～12.3%.粉碎區(qū)半徑隨著埋深的變化幅度較小，其主要是爆炸產(chǎn)生的能量主要向自由面方向傳播，而埋深對爆轟波以及爆生氣體的膨脹作用影響較小所導(dǎo)致.

碎裂區(qū)形狀呈類水滴狀，炮孔軸向的碎裂區(qū)深度在757～941 mm 之間，超出炮孔深度34.4%～51.4%，巖石以壓縮破壞為主，基本無自由面效應(yīng)；炮孔徑向的破碎區(qū)半徑隨著埋深的增加存在一定規(guī)律，在數(shù)值模擬結(jié)果中提取距自由面不同深度破碎區(qū)范圍見圖7，可知：接近自由面位置碎裂區(qū)半徑值最大，之后隨著深度增加，呈現(xiàn)下降趨勢；通過結(jié)果數(shù)據(jù)擬合后獲得了碎裂區(qū)半徑rds隨自由面距離df的變化規(guī)律為

圖7 碎裂區(qū)半徑變化趨勢Fig.7 Radius variation trend of fracture zone

（3）爆破漏斗體積變化規(guī)律

對各組別的爆破漏斗進行可視化處理并計算其體積，結(jié)果如圖8 所示，可以看出：爆破漏斗體積呈先增大后減小的趨勢，在埋深為590 mm 時達到最大.數(shù)據(jù)顯示，粉碎區(qū)半徑為裝藥半徑的18.8～30.6倍，碎裂區(qū)半徑為裝藥半徑的36.2～40.4倍，損傷區(qū)半徑為裝藥半徑的51.3～70 倍，符合爆破動力學(xué)[46]中巖石爆破的碎裂區(qū)、粉碎區(qū)半徑確定公式的理論值.

圖8 爆破漏斗模擬體積Fig.8 Simulation volume of blasting funnel

3 碎裂區(qū)范圍現(xiàn)場爆破試驗驗證

3.1 現(xiàn)場爆破試驗情況

為驗證數(shù)值模擬過程中RHT 本構(gòu)碎裂區(qū)判定范圍的可靠性，在大紅山鐵礦400 m 中段8 號穿脈巷道側(cè)壁進行多組變埋深爆破漏斗試驗，鉆鑿工作采用氣腿式鑿巖機，有效控制偏斜與鉆孔深度誤差，設(shè)計試驗按照數(shù)值仿真方案中的參數(shù)，共開展5 組，每組3 個孔，孔間間隔2.5 m 以上，爆破漏斗試驗場地以及所用爆破器材如圖9 所示.

圖9 爆破漏斗試驗現(xiàn)場及所用設(shè)備Fig.9 Blasting funnel test site and equipment

3.2 爆破結(jié)果分析

爆破后采用先進的高精度Maptek SR3 三維激光掃描儀對爆破漏斗進行掃描，通過后處理軟件對爆破漏斗體積進行可視化處理，爆破漏斗試驗效果及數(shù)據(jù)統(tǒng)計分析結(jié)果如圖10 所示.由試驗結(jié)果可知：爆破漏斗形狀較不規(guī)律，部分存在大塊脫落的現(xiàn)象，爆破漏斗半徑約為裝藥半徑的22～30 倍；同數(shù)值模擬結(jié)果相比存在普遍偏大現(xiàn)象，各項平均擬合度分別為91.3%、89%、81.3%，整體擬合度較高，證明RHT碎裂區(qū)判定范圍具有一定的可靠性.

圖10 爆破漏斗數(shù)據(jù)對比分析Fig.10 Comparative analysis of blasting funnel data

3.3 破碎區(qū)優(yōu)化

將圖中數(shù)據(jù)進行對比可以看出，數(shù)值模擬所劃分的半徑范圍與現(xiàn)實所測半徑相差較大，因此對式（11）進行細分優(yōu)化，優(yōu)化后碎裂區(qū)半徑rd、粉碎區(qū)半徑rc與損傷變量D之間的關(guān)系為

優(yōu)化后的碎裂區(qū)、粉碎區(qū)半徑確定理論公式能夠更好地為大紅山鐵礦相關(guān)爆破工程參數(shù)設(shè)計提供理論計算依據(jù).

4 結(jié) 論

針對RHT 本構(gòu)模型在損傷判定相關(guān)研究方面的部分缺陷，本文基于巖石損傷準則與RHT 本構(gòu)方程建立巖石爆破損傷分區(qū)的判定方法，基于物質(zhì)點法對變埋深條件下的爆破漏斗成型規(guī)律開展系列數(shù)值仿真，并結(jié)合多組現(xiàn)場試驗對巖石碎裂區(qū)、粉碎區(qū)進行驗證，得出以下主要結(jié)論：

（1）基于巖石損傷判定準則與RHT 內(nèi)嵌本構(gòu)方程對巖石爆破過程造成的損傷分區(qū)進行理論推導(dǎo)，計算出以變鈉質(zhì)熔巖為主要巖性下的損傷分區(qū)臨界閾值Dcr=0.11，碎裂分區(qū)臨界閾值Dcf=0.51，并基于損傷半徑確定公式求出不同損傷區(qū)爆心距.

（2）利用物質(zhì)點法對變埋深爆破漏斗試驗進行數(shù)值計算，獲得了變鈉質(zhì)熔巖爆破漏斗成型過程中不同時間不同區(qū)域最大應(yīng)變量；不同埋深條件下，粉碎區(qū)、破碎區(qū)及損傷區(qū)范圍和半徑；碎裂區(qū)半徑隨自由面距離的變化規(guī)律.

（3）現(xiàn)場爆破漏斗試驗驗證了RHT 本構(gòu)碎裂區(qū)判定范圍具有一定的可行性；對碎裂區(qū)、粉碎區(qū)損傷半徑確定公式進行了優(yōu)化，獲得了變鈉質(zhì)熔巖碎裂區(qū)、粉碎區(qū)半徑與損傷變量之間的關(guān)聯(lián)關(guān)系.