朱心廣，馮 春，王心泉，程鵬達(dá)，高圣元

（1. 中國科學(xué)院力學(xué)研究所，北京 100190；2. 中國科學(xué)院大學(xué)工程科學(xué)學(xué)院，北京 100049；3. 煤炭科學(xué)研究總院，北京 100013）

臺階爆破技術(shù)作為一種高效經(jīng)濟(jì)的破巖手段，廣泛應(yīng)用于交通、礦山、水利水電等各個(gè)領(lǐng)域[1]。隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展，數(shù)值分析在臺階爆破技術(shù)中參數(shù)的選擇及優(yōu)化方面扮演著越來越重要的角色。目前關(guān)于爆炸問題的數(shù)值計(jì)算方法，按照有無網(wǎng)格可分為兩類。網(wǎng)格類計(jì)算方法包括有限元法[2-5]和離散元法[6-8]，其中有限元法的數(shù)學(xué)模型和算法相對成熟，針對連續(xù)問題的分析以及應(yīng)力場的描述方面具有優(yōu)勢，但是在大位移、大變形問題中常因單元過分畸變導(dǎo)致計(jì)算發(fā)散；離散元法在描述塊體的大位移、碰撞、飛散等過程中具有優(yōu)勢，但是無法描述連續(xù)損傷過程。在網(wǎng)格類的計(jì)算方法中，為了獲得精確結(jié)果，需要在爆源附近加密網(wǎng)格，導(dǎo)致網(wǎng)格數(shù)量增大，增加了計(jì)算時(shí)間。無網(wǎng)格的計(jì)算方法包括物質(zhì)點(diǎn)法[9-10]（material point method, MPM）、光滑粒子流體動力學(xué)[11-14]（smoothed particle hydrodynamics, SPH）等，其突出特點(diǎn)是在計(jì)算過程中與網(wǎng)格無關(guān)，因此不存在網(wǎng)格加密的問題，但是和有限元相比數(shù)學(xué)模型和算法還不夠成熟，計(jì)算精度和計(jì)算效率還有一定差距。

對于爆破問題的數(shù)值模擬，比較成熟的商業(yè)軟件包括LS-Dyna、Autodyn、UDEC 和3DEC 等，常用的計(jì)算方法包括非連續(xù)變形分析法（discontinuous deformation analysis, DDA）、有限離散元（finite discrete element method, FDEM）以及連續(xù)-非連續(xù)單元法（continuous-discontinuous element method, CDEM）等。余仁兵等[15]通過LS-Dyna 軟件對不同直徑炮孔不耦合裝藥結(jié)構(gòu)簡化模型進(jìn)行了數(shù)值模擬；丁希平[16]采用LS-Dyna 軟件研究了炸藥單耗、超深及填塞長度對爆破效果的影響；璩世杰等[17]運(yùn)用LS-Dyna 軟件研究了不同角度的節(jié)理對預(yù)裂爆破效果的影響。劉超[18]使用Autodyn 軟件針對兩種不同的間隔裝藥方式建立計(jì)算模型，計(jì)算結(jié)果表明中間隔位置錯(cuò)開的裝藥結(jié)構(gòu)相對中間隔位置對齊的裝藥結(jié)構(gòu)而言，各處壓力大小分布更加均勻，有效作用時(shí)間更長；謝冰等[19]運(yùn)用Autodyn 與UDEC 相結(jié)合的方法，研究了節(jié)理幾何特征對預(yù)裂爆破的影響。Yan 等[20]利用3DEC 軟件，通過在炮孔破碎圈外側(cè)施加等效破碎載荷，并控制不同區(qū)域的網(wǎng)格尺寸，研究了三維條件下爆堆的形成過程。郭雙等[21]基于DDA 方法建立了巖石爆破力學(xué)模擬模型，分別研究了地應(yīng)力對爆生氣體壓力和爆炸應(yīng)力波破巖效果的影響；江成等[22]采用DDA 方法模擬了爆炸應(yīng)力波作用下考慮地應(yīng)力條件時(shí)的單孔和多孔鑿巖爆破破巖過程。Han 等[23]采用FDEM 方法并運(yùn)用GPU 并行加速的方式，研究了爆破作用下，對某含有臺階段深埋隧道的損傷演化過程。鄭炳旭等[24]基于CDEM 方法，探討了炸藥單耗對爆破塊度分布曲線及系統(tǒng)破裂度的影響規(guī)律；馮春等[25]通過在CDEM 方法中引入朗道爆源模型、彈性-損傷-斷裂本構(gòu)及半彈簧接觸算法[26]，模擬了從炸藥起爆、巖體損傷破裂到最后爆堆形成的全過程。

在爆源模型的處理上，ABAQUS 和LS-Dyna 通常建立炸藥與空氣數(shù)值分析模型并借助其狀態(tài)方程施加爆炸荷載，該方法需要對炸藥和空氣分別單獨(dú)劃分網(wǎng)格并賦予其各自的狀態(tài)方程，且炸藥網(wǎng)格和空氣網(wǎng)格應(yīng)在交界面上共節(jié)點(diǎn)，操作較為復(fù)雜且需要較長的計(jì)算時(shí)間。一般而言，鉆孔爆破的炮孔尺寸在厘米至分米量級，而待爆巖體的范圍則在米至數(shù)十米范圍。數(shù)值模擬時(shí)，若將炮孔與待爆巖體均用實(shí)體網(wǎng)格進(jìn)行剖分，將會導(dǎo)致炮孔附近的網(wǎng)格過密，從而使得整個(gè)模型的網(wǎng)格數(shù)量過大而影響計(jì)算效率。為解決上述問題，本文基于CDEM[27]方法，提出一種利用線狀桿件表述炮孔及炸藥，通過線狀桿件與周圍實(shí)體單元的耦合實(shí)現(xiàn)爆破破巖過程模擬的計(jì)算方法。

1 桿件爆源模型

1.1 基本假設(shè)

(1) 炸藥爆轟過程中滿足軸對稱假設(shè)，爆炸過程中炸藥桿件單元在軸向上不發(fā)生膨脹，僅在徑向發(fā)生膨脹。

(2) 炸藥與實(shí)體單元滿足壓力-體積耦合關(guān)系，即炸藥單元傳遞給實(shí)體單元膨脹應(yīng)力，實(shí)體單元傳遞給炸藥膨脹體積。

(3) 不考慮爆生介質(zhì)與炮孔壁對的相互作用，直接將爆轟產(chǎn)物產(chǎn)生的氣體壓力作用于周邊圍巖，以加快求解速度。

由于不考慮爆轟產(chǎn)物與炮孔的相互作用過程，因此該模型對炸藥近場的模擬不是很準(zhǔn)確。如圖1 所示，Ei表示待爆巖體網(wǎng)格，Ni表示桿件單元節(jié)點(diǎn)（i=1,2,3,···）。實(shí)施本文所述的方法時(shí)，首先將待爆巖體采用較大的實(shí)體網(wǎng)格進(jìn)行剖分，然后將炮孔簡化為桿件插入至待爆巖體的指定位置；桿件上劃分若干單元，單元中引入經(jīng)典的朗道爆源模型；判斷每個(gè)桿件節(jié)點(diǎn)與實(shí)體單元的拓?fù)潢P(guān)系，若該桿件節(jié)點(diǎn)位于某個(gè)實(shí)體單元的內(nèi)部或表面，則將該實(shí)體單元作為桿件節(jié)點(diǎn)的傳力對象，將桿件節(jié)點(diǎn)上的爆生氣體壓力施加至該巖石單元上；同時(shí)，根據(jù)該實(shí)體單元的體應(yīng)變計(jì)算出桿件節(jié)點(diǎn)處的橫截面變化情況，用于計(jì)算下一時(shí)刻的爆生氣體壓力。

圖1 桿件單元與實(shí)體網(wǎng)格的耦合示意Fig. 1 Schematic of coupling of bar and solid mesh

1.2 朗道爆源模型

朗道模型采用朗道-斯坦紐科維奇公式（γ 率方程）進(jìn)行爆炸氣體膨脹壓力的計(jì)算：

式中：γ1及γ2分別為第一段及第二段的絕熱指數(shù)，對于一般的凝聚態(tài)炸藥，γ1= 3，γ2= 4/3；p和V分別為高壓氣球的瞬態(tài)壓力和體積，p0和V0分別為高壓氣球初始時(shí)刻的壓力和藥包的體積，pc和Vc分別為高壓氣球在兩段絕熱過程邊界上的壓力和體積。pc和p0分別為：

式中：Qw為炸藥爆熱，J/kg；ρw為裝藥密度，kg/m3；D為爆轟速度，m/s。

此外，當(dāng)某時(shí)刻爆炸產(chǎn)生的壓力大于CJ (Chapman-Jouguet) 面上的壓力 (pCJ) 時(shí)，將該壓力修正為CJ 面上的壓力。

本文采用線狀桿件描述炸藥，每個(gè)炸藥單元的初始體積及當(dāng)前體積為

式中：r0為炸藥單元的初始半徑，h為炸藥單元的長度，θe為炸藥單元的平均體應(yīng)變 (可通過與該單元存在耦合關(guān)系的巖石單元的體應(yīng)變獲得)。

為了模擬點(diǎn)火爆炸及爆轟波的傳播過程，朗道模型中需設(shè)置點(diǎn)火點(diǎn)位置、點(diǎn)火時(shí)間，并根據(jù)到時(shí)起爆判斷某一炸藥單元是否執(zhí)行爆炸壓力計(jì)算。

設(shè)某一炸藥 (含若干個(gè)桿件單元) 的點(diǎn)火時(shí)間為t0，點(diǎn)火點(diǎn)坐標(biāo)為(x0,y0,z0)，該炸藥中某一桿件單元的中心到點(diǎn)火點(diǎn)的距離為d，炸藥的爆速為D，則該單元的點(diǎn)火時(shí)間為t1=d/D+t0。當(dāng)爆炸時(shí)間t＞t1時(shí)，該單元開始進(jìn)行爆炸壓力的計(jì)算，計(jì)算公式為

式中：pr為真實(shí)爆炸壓力，f(p)為爆轟產(chǎn)物狀態(tài)方程 (可根據(jù)式(1)獲得); ζ 為能量釋放率：

鉆孔爆破過程中，隨著炮孔的起爆，爆生氣體逐漸膨脹，并推動巖體做功。巖體在爆炸壓力的作用下逐漸出現(xiàn)裂縫的萌生、擴(kuò)展，并最終出現(xiàn)貫通裂縫。一旦出現(xiàn)貫通性的裂縫，炮孔內(nèi)的氣體將會從裂縫中快速溢出，并導(dǎo)致炮孔內(nèi)的壓力急劇下降。由于直接模擬爆生氣體在巖體內(nèi)的流動及溢出過程較為復(fù)雜，本文采用設(shè)置失效時(shí)間的方式進(jìn)行等效模擬。當(dāng)爆破時(shí)間大于失效時(shí)間時(shí)，該炮孔即設(shè)置為失效單元，失效單元中的氣體壓力pr為0，且不再進(jìn)行爆炸壓力的計(jì)算。

1.3 桿件單元與巖石單元的耦合

設(shè)某線狀炸藥中共有N個(gè)節(jié)點(diǎn)位于某個(gè)實(shí)體單元內(nèi)部或表面 (圖1)，則該線狀炸藥位于實(shí)體單元內(nèi)部的平均爆壓可表述為

式中：負(fù)號是因?yàn)楸▔毫σ詨簽檎?，而?yīng)力以拉為正；σb為施加于實(shí)體單元上的膨脹應(yīng)力，即可認(rèn)為是施加在實(shí)體單元上的三項(xiàng)正應(yīng)力（剪應(yīng)力為0）；α為壓力衰減指數(shù)（需根據(jù)數(shù)值實(shí)驗(yàn)進(jìn)行標(biāo)定）；d為實(shí)體單元的特征尺寸：

式中：Ae為二維實(shí)體單元的面積，Ve為三維實(shí)體單元的體積。

將上述膨脹應(yīng)力轉(zhuǎn)化為節(jié)點(diǎn)力施加于單元的各個(gè)節(jié)點(diǎn)上，為

式中：fni為單元中節(jié)點(diǎn)n第i個(gè)方向的節(jié)點(diǎn)力分量，M為單元中與節(jié)點(diǎn)n相關(guān)的面的數(shù)量，eki為與第k個(gè)面的單位外法向量在第i個(gè)方向的分量，Ak為第k個(gè)面的面積，Qk為第k個(gè)面的節(jié)點(diǎn)總數(shù)。

式(4)中炸藥單元的平均體應(yīng)變θe可表示為

式中：p、q為炸藥單元的兩個(gè)節(jié)點(diǎn)，θp及θq為與節(jié)點(diǎn)p、q 存在耦合關(guān)系的實(shí)體單元的體應(yīng)變。

因此，雖然假定了炸藥單元僅在徑向發(fā)生膨脹，但是通過膨脹壓力以正應(yīng)力形式施加于實(shí)體單元節(jié)點(diǎn)，也可以實(shí)現(xiàn)實(shí)體單元的軸向膨脹破壞。

2 壓力衰減指數(shù)α 的標(biāo)定

在半無限巖體內(nèi)設(shè)置一個(gè)直徑為250 mm 的鉆孔，鉆孔深度為15 m，鉆孔內(nèi)的裝藥長度為8 m，填塞7 m。分別建立實(shí)體單元數(shù)值模型(模型A)及桿件-實(shí)體單元混合數(shù)值模型(模型B)，通過對比起爆后典型位置的徑向幅值振動速度vp，來標(biāo)定式(9)中壓力衰減指數(shù)α 的具體數(shù)值。

由于在實(shí)體炮孔模型中，要求實(shí)體炮孔和圍巖共節(jié)點(diǎn)，而桿件模型中不要求桿件與圍巖單元共節(jié)點(diǎn)，從這一點(diǎn)來講，無法做到兩種數(shù)值模型完全一致。另外，本節(jié)的目的就是改變不同網(wǎng)格尺寸標(biāo)定α 值，只需保證測點(diǎn)與爆心的距離相等即可，因此選擇什么樣的幾何構(gòu)型不是很重要。

如圖2 所示，模型A 為1/4 圓柱模型，高為25 m，半徑為25 m，炸藥位于柱體中心軸線上。模型的外周面及底部設(shè)置為無反射邊界條件，模型的對稱面上施加法向位移約束條件。采用映射網(wǎng)格進(jìn)行剖分，軸向共分割50 份 (每份0.5 m)，環(huán)向共分割15 份，炮孔部分的徑向共分割3 份，巖體中的徑向采用8、16、48 三種分割數(shù)進(jìn)行分割，所建立的三種網(wǎng)格分別命名為網(wǎng)格A1 (每份3.11 m)、網(wǎng)格A2 (每份1.55 m)、網(wǎng)格A3 (每份0.518 m)。

模型B 的炸藥采用桿件單元進(jìn)行剖分，巖體采用六面體單元進(jìn)行剖分 (如圖3 所示)。巖體的高度25 m，長度及寬度均為50 m，巖體高度方向共分割25 份 (每份1 m)，水平方向分別采用16、25、32、50 四種分割數(shù)進(jìn)行分割，所建立的4 種網(wǎng)格分別命名為網(wǎng)格B1 (每份3.125 m)、網(wǎng)格B2 (每份2 m)、網(wǎng)格B3 (每份1.5625 m)、網(wǎng)格B4 (每份1 m)。炸藥單元位于模型正中，采用42 個(gè)桿件單元進(jìn)行離散。為了消除人工邊界的虛假反射，在模型四周及底面均施加無反射邊界條件。

本節(jié)中巖體為赤鐵礦，采用Mohr-Coulomb理想彈塑性模型(巖體單元間不考慮斷裂)進(jìn)行描述；炸藥為現(xiàn)場混裝乳化炸藥，采用朗道爆源模型進(jìn)行描述。巖體參數(shù)及炸藥參數(shù)見表1 及表2。

表1 鐵礦石的力學(xué)參數(shù)Table 1 Mechanical parameters of iron ore

表2 乳化炸藥的參數(shù)Table 2 Parameters of emulsion

為了標(biāo)定壓力衰減指數(shù) α ，選取模型A 及模型B 中的對應(yīng)測點(diǎn)進(jìn)行對比(測點(diǎn)位置見圖2、圖3)。該測點(diǎn)到爆源的徑向距離為25 m，兩個(gè)模型中測點(diǎn)的徑向振動速度峰值隨單元尺寸的變化如圖4 所示。由圖4 可得，隨著單元尺寸的增大，兩個(gè)模型所獲得的測點(diǎn)處的徑向速度幅值均線性減??；當(dāng)壓力衰減指數(shù) α 為1.25 時(shí)，模型B 所獲得的vp衰減規(guī)律與模型A 所獲得的規(guī)律基本一致。由此表明，當(dāng)式(8)中的壓力衰減指數(shù)取1.25 時(shí)，一維軸對稱爆源模型的爆破效果與實(shí)體模型的爆破效果基本一致。

圖2 模型A (徑向分割16 份)Fig. 2 Model A (dividing into 16 segment in the radial direction)

圖3 模型B (水平向分割16 份)Fig. 3 Model B (dividing into 16 segment in the horizontal direction)

圖4 測點(diǎn)處徑向峰值振動速度隨單元尺寸的變化Fig. 4 Variation of radial peak particle velocity at monitoring points with element size

網(wǎng)格A2 及網(wǎng)格B3 對應(yīng)的徑向網(wǎng)格尺寸均為1.5 m 左右，當(dāng)模型B 中軸對稱爆源模型所用的壓力衰減指數(shù) α 為1.25 時(shí)，兩個(gè)網(wǎng)格獲得的測點(diǎn)處的徑向振動速度時(shí)程曲線如圖5 所示。由圖5 可得，兩者的振動規(guī)律基本一致，僅起振時(shí)間及首峰后的振動規(guī)律有所差別，由此進(jìn)一步證明了衰減指數(shù) α 為1.25 時(shí)的一維軸對稱爆源模型，其爆破效果可以與實(shí)體模型的效果基本一致。

圖5 測點(diǎn)處徑向振動速度時(shí)程曲線Fig. 5 History of radial velocity at monitoring point

3 模型驗(yàn)證

建立如圖6 所示的立方體混凝土模型，內(nèi)部設(shè)置炮孔模型，混凝土尺寸為3.048 m×3.048 m×1.52 m，炮孔直徑38 mm，長1200 mm，炮孔位置見剖面圖?；炷良芭诳琢W(xué)參數(shù)采用文獻(xiàn)[28]中數(shù)據(jù)，其中混凝土密度為2009 kg/m3，抗拉強(qiáng)度為1.8 MPa，彈性模量為13.1 GPa，泊松比為0.25；炮孔直徑為38 mm，密度1150 kg/m3，爆速4700 m/s1，爆熱3.24×106J/kg，桿件爆源中爆炸持續(xù)時(shí)間為15 ms。

當(dāng)模型爆炸20 ms 時(shí)，位于炮孔中部位置的模型爆炸結(jié)果剖視圖如圖7 所示，其中最大位移發(fā)生在炮孔左側(cè)的自由臨空面處，最大位移值0.88 m，整個(gè)爆炸形態(tài)與文獻(xiàn)[28]結(jié)果基本一致。

圖7 位于炮孔中部位置的模型爆炸結(jié)果剖視圖(20 ms)Fig. 7 Sectional view of model explosion results in the middle of the blast hole after 20 ms

4 工程應(yīng)用

鞍鋼露天鐵礦臺階爆破現(xiàn)場圖片如圖8 所示。本文以鞍鋼露天鐵礦臺階爆破為例，根據(jù)爆破設(shè)計(jì)（圖9(a)），建立如圖9(b)所示的5 排10 列 (共50 個(gè)炮孔) 數(shù)值模型，該模型共包含19658 個(gè)計(jì)算網(wǎng)格，其中四面體單元數(shù)6541，三棱柱單元數(shù)720，金字塔單元數(shù)6916，六面體單元數(shù)6381. 根據(jù)現(xiàn)場情況，模型表面包含3 m 厚的松散堆積石土，下部為鐵礦石塊。炮孔深度為15 m，直徑為0.25 m，超深3 m，段高12 m，間排距為7 m×7 m，單孔裝藥量677 kg，采用逐孔起爆方式，孔底起爆，孔間延時(shí)42 ms，排間延時(shí)65 ms，連線方式為串并聯(lián)聯(lián)結(jié)。

圖8 鞍鋼露天礦現(xiàn)場Fig. 8 Site view of Angang strip mine

圖9 含50 炮孔的臺階爆破數(shù)值模型Fig. 9 The bench blasting model with 50 bore holes

松散堆積石土及鐵礦石塊采用Mohr-Coulomb 理想彈塑性模型進(jìn)行描述，其中松散堆積石土的密度為1800 kg/m3，彈性模量為6 GPa，泊松比為0.25，內(nèi)聚力15 MPa，抗拉強(qiáng)度5 MPa，內(nèi)摩擦角40°；鐵礦石塊彈性模量為4 GPa，其他力學(xué)參數(shù)見表1。爆炸所用炸藥采用現(xiàn)場混裝的乳化炸藥，其參數(shù)見表2。

數(shù)值計(jì)算共分為兩個(gè)階段：第一階段采用虛擬質(zhì)量法獲得模型在重力作用下的靜態(tài)應(yīng)力場，該階段中模型四周及底部為法向約束邊界，重力方向豎直向下；第二階段起爆階段，模型底部及四周施加無反射邊界，計(jì)算時(shí)步為0.3 μs，Rayleigh 阻尼中剛度阻尼系數(shù)為1×10-4，質(zhì)量阻尼系數(shù)為0。

臺階爆破后露天邊坡的破壞狀態(tài)如圖10 所示，由圖可以直觀看出，爆區(qū)內(nèi)整體以拉伸破壞為主，只有在爆區(qū)頂部及臺階拐角處出現(xiàn)零星剪切破壞。

圖10 臺階爆破后邊坡的破壞狀態(tài)Fig. 10 Failure status of slope after bench blasting

在地表布設(shè)5 個(gè)監(jiān)測點(diǎn)，分別為M1～M5(圖9(b))，其中M1 與M2 的間距為18.5 m，M2 與M3 的間距為17.2 m，M3 與M4 的間距為17.8 m，M4 與M5 的間距為18.5 m。獲取各個(gè)測點(diǎn)處振動速度峰值并與現(xiàn)場實(shí)測數(shù)據(jù)對比，繪制曲線如圖11 所示。由圖11 可得，爆點(diǎn)與測點(diǎn)的水平距離越大，測點(diǎn)的振動峰值速度越小，數(shù)值計(jì)算結(jié)果與現(xiàn)場實(shí)測數(shù)據(jù)的誤差越小，說明運(yùn)用本文所提出的桿件爆源模型進(jìn)行遠(yuǎn)場的爆炸模擬的結(jié)果是可信的，也證明了桿件爆源模型的有效性與合理性。

圖11 峰值振動速度隨距離的變化曲線Fig. 11 The curve of peak particle velocity with distance

5 結(jié) 論

(1) 提出了一種一維軸對稱桿件爆源模型，利用線狀桿件表述炮孔及炸藥，炸藥(桿件)單元傳遞給周圍實(shí)體單元膨脹應(yīng)力，實(shí)體單元傳遞給炸藥膨脹體積，兩者通過壓力-體積耦合關(guān)系實(shí)現(xiàn)炸藥與周圍巖體的相互作用。

(2) 通過與實(shí)體炮孔模型的數(shù)值對比分析，當(dāng)壓力衰減指數(shù) α =1.25 時(shí)，本文所提爆源模型獲得的徑向振動速度峰值衰減規(guī)律與實(shí)體炮孔模型所獲得的規(guī)律基本一致，并且針對混凝土塊破壞特性分析，其破壞形態(tài)與文獻(xiàn)基本一致。

(3) 以鞍鋼露天鐵礦臺階爆破開采為背景，模擬了5 排10 列共50 個(gè)炮孔逐孔起爆后，爆區(qū)內(nèi)的損傷破壞狀態(tài)。數(shù)值計(jì)算結(jié)果表明，爆區(qū)內(nèi)以拉伸破壞為主，并且所獲得測點(diǎn)處的振動速度峰值隨距離的變化規(guī)律與現(xiàn)場實(shí)測數(shù)據(jù)基本一致。

(4) 采用一維軸對稱爆源模型后，炮孔附近不再需要加密網(wǎng)格，因此可以降低計(jì)算網(wǎng)格數(shù)量并增大計(jì)算步長，并最終實(shí)現(xiàn)炸藥起爆過程及爆炸地震波傳播過程的統(tǒng)一高效模擬。