土中爆炸成腔在地下工程、军事爆破、地震勘探等方面有重要的应用前景。一般地震勘探中使用的炸药震源因为药量大、埋藏深的原因，球形药包不易装药，工程上往往使用柱形装药的方式，其长径比为10～20，不同的装药方式会产生不同的爆腔，进而影响远场地震波的频谱特征。

但是由于土介质参数及土壤结构的复杂性，其爆腔范围不易准确计算。谢多夫^[1]提出了理想介质爆腔半径的计算方法，但该方法由理想空气介质中的爆炸过程推导而来，在土介质中使用误差较大。王海亮等^[2]、张奇^[3]对柱形装药爆腔膨胀过程进行了研究，给出了爆腔与装药半径之间的近似倍数关系。王仲琦等^[4]对爆炸挤压黏土产生的密度变化过程进行了数值模拟。目前，针对土中爆炸成腔半径的理论计算方法有基于爆炸相似理论的半经验公式^[5]、准静定理论爆炸成腔计算方法^[5]、条形装药爆破成腔半径的弹塑性估算^[6]、球形装药爆腔准静态模型^[7]、柱形装药爆腔准静态预测模型^[8]等。但是，半经验公式和弹塑性估算计算误差较大；准静定理论所需参数较多、计算复杂；柱形装药爆腔准静态模型只能计算径向尺寸，不能计算轴向尺寸。能够快速给出柱形药包爆腔特征尺寸的计算方法较少。

在球形装药爆腔预测准静态模型提出之前，柱形药包爆腔理论计算参数复杂、适用范围有限，而且利用球形药包代替作近似计算也较困难^[9]。本文中，以球形装药爆腔预测准静态模型为基础，对长径比较大的柱形药包进行近似替代，并给出柱形药包爆腔和塑性区的特征尺寸，和数值模拟结果进行比较分析，以期研究结果对解决地震勘探工程中柱形药包的成腔问题起到参考作用。

1.   柱形爆腔近似计算方法的建立

1.1   柱形和球形装药爆腔预测准静态模型

土介质相对松软，在强荷载作用下发生剧烈变形，由于炸药爆炸产生爆生气体发生膨胀，当爆炸冲击波波阵面压力远高于土体强度时，土介质颗粒与冲击波同向位移，水和气体被挤出形成爆腔。在爆轰过程中，炸药产生高温高压气体产物并对周围土介质作用，初始阶段爆生气体产物压力较高，产生的应力波波阵面压力远高于土的屈服极限，土体呈流体性质发生剧烈压缩变形形成空腔（空腔区），当应力波进入空腔外的土介质中，土体受挤压吸收能量，应力低于土体破坏极限但高于弹性极限，发生塑性变形（塑性区），随着应力波在塑性区的吸收、扩散，衰减，应力低于土体弹性极限，土体发生弹性变形（弹性区），应力波变为弹性波。

Drukovanyi等^[8]对土中爆炸特征进行了丰富的研究，并给出了一维柱形药包爆腔准静态理论模型。假设柱形药包和周围土介质在轴向和径向都无限延伸，空腔膨胀过程是轴对称的，切向应力为零，剪切应变和轴向应变也为零，引入黏聚力颗粒状不可压缩介质模型^[9]，得出爆腔半径和塑性区半径。

爆腔半径为：

塑性区半径为：

式中：a₀为球形药包半径，p₀为炸药初始爆压，γ为炸药膨胀指数，f为土介质内摩擦角，c为土介质黏聚力，σ_*为土介质抗压强度，σ₀为土介质抗拉强度，μ为拉梅系数(在弹性介质中为剪切模量)。

于成龙等^[7]基于一维柱形装药准静态模型建立了球形装药爆腔预测准静态模型，并作出3种假设条件：(1)岩土介质不可压缩；(2)炸药爆炸对岩土介质的改变忽略不计；(3)炸药爆轰过程和空腔的形成都是瞬时。同时推出球形装药爆腔半径和塑性区半径。

爆腔半径为：

塑性区半径为：

1.2   计算方法

根据上述球形装药准静态模型，利用球形叠加的方式，提出一种计算柱形装药爆腔近似方法(假设条件一致)，当柱形装药长径比足够大时(本文中取长径比大于等于10)，可近似认为该柱形药包由n个合适的球形药包叠加组成。在炸药性质和质量不变的情况下，应满足：

（1）柱形药包和球形药包总体积相等；

（2）爆轰过程的连续性，球形药包间距为零；

（3）药包形状相近，球形药包直径应尽量与柱形药包直径相近或相等，球形药包直径之和应尽量与柱形药包长度相近或相等。那么可分为3种近似替代情况，如图1所示。

图  1  划分方式

Figure  1.  Division methods

下载: 全尺寸图片 幻灯片

为便于求解，优先考虑球形药包直径之和与柱形药包长度相等的情况，即图1(a)的情况：

式中：r₁、r₂分别为柱形药包半径和球形药包半径，且2r₁＞r₂＞r₁。

同理，图1(b)的解为：

图1(c)的解为：

则柱形药包爆腔的特征尺寸可用以下公式表示。

爆腔长半轴为：

爆腔短半轴仍为b_*。

塑性区长半轴为：

塑性区短半轴仍为b₀，其中N代表球形药包个数n、n+1、n+2等3种情况中的某一值。

为比较3种近似替代情况的合理性，对6 kg TNT在土中爆炸形成的爆腔进行计算分析。土介质参数^[10]为：抗压强度σ_*=13 MPa，抗拉强度σ₀=2 MPa，拉梅系数μ=0.147 GPa，内摩擦角f=0.115，黏聚力c=11.8 kPa，密度ρ_soil=1 840 kg/m³，体积模量K=245 MPa，剪切模量G=147 MPa，屈服强度σ_s=22 MPa。

根据式(6)～(12)计算得出图1(a)、(b)、(c)分别对应的球形药包个数N值为8、9、10，对应的球形药包半径分别为47.89、46.05、44.46 mm。根据式(1)、(13)球形药包近似替代柱形药包与柱形药包形成的爆腔^[8]计算结果如图2所示。

图  2  不同N值的爆腔包络图(单位：mm)

Figure  2.  Burst envelope diagrams at different N values (unit: mm)

下载: 全尺寸图片 幻灯片

从图2的爆腔包络图中可以看出，柱形药包转球形药包划分的球形药包个数越多，球形药包的直径之和越大，形成的爆腔轴向越长，包络图越扁长，与一维柱形装药准静态模型半径相差越大。而N=8时，两者爆腔半径最接近，即按照式(13)计算当N=n时，图1(a)的近似替代方法较合理，此时${b_*} = 189.5$ mm，${b_{*n}} = 524.73$ mm。

2.   数值模拟验证分析

为比较3种近似替代方法的合理性，使用显式有限元分析程序AUTODYN进行二维数值模拟。AUTODYN可用来解决固体、流体、气体及其相互作用的高度非线性动力学问题，广泛应用于冲击、爆炸等问题的分析研究。本文中，用该软件模拟土介质材料在不同形状药包作用下的成腔过程。

2.1   炸药爆轰产物状态方程及参数

炸药选择TNT，采用标准JWL理想爆轰产物状态方程，其表达式如下：

式中：e为CJ爆轰能量，C₁、C₂、R₁、R₂、ω为试验确定常数，V为相对体积。试验确定常数C₁=373.7 GPa，C₂=3.747 GPa，R₁=4.15，R₂=0.90，ω=0.35；炸药密度ρ=1 650 kg/m³，爆轰速度D=6 930 m/s，CJ爆轰能量e=6.0 GJ/m³，CJ压力p_CJ=21 GPa，初始爆压p₀=9.82 GPa。

2.2   土介质物理模型

土介质采用理想弹塑性模型，状态方程为：

式中：K为土介质的体积模量。

为体现土介质在爆炸作用下的弹塑性变化过程，土介质的强度模型则采用von Mises强度模型：

式中：σ_s为土介质的屈服强度，G为介质的剪切模量。

土介质使用1.2节中所给出的参数。

2.3   几何模型

取6 kg TNT炸药，柱形药包长径比为10，按照表1的2种方式划分成柱形和球形药包进行模拟。对于炸药及炸药近区，因网格变形过大，使用Euler单元网格描述；远区冲击波应力低于介质强度极限，变形较小，则使用Lagrange单元描述；接触部分使用Euler/Lagrange耦合算法；土介质远区使用transmit完全投射边界条件，左右对称轴处使用velocity条件，x方向速度为零。因爆腔最终结果与药包起爆点位置无关^[11]，爆腔形状存在对称性，因此使用四分之一网格进行计算，减少计算量。Euler单元网格介质尺寸为1 100 mm×500 mm，网格划分为110×50，Lagrange单元网格介质尺寸为4 000 mm×4 000 mm，网格划分为200×200。

根据参考文献[8]，柱形药包爆腔近似为椭球型，在二维条件下，用离心率表征比较两种划分方式的近似程度，计算得出爆腔特征尺寸。计算结果如表1和图3所示。

表  1  两种药包形成的爆腔特征尺寸

Table  1.  Characteristic dimensions of blasting zones formed by two kinds of charge

划分方式	等体积球 个数	半径/ mm	爆腔长半轴/ mm	爆腔短半轴/ mm	爆腔离心率	塑性区长半轴/ mm	塑性区短半轴/ mm	塑性区离心率
柱形		38.84	530	200	0.377 4	800	570	0.712 5
球形	8	47.89	500	190	0.380 0	790	560	0.708 9

下载: 导出CSV 

| 显示表格

图  3  两种药包形成的爆腔特征尺寸图(深蓝色为爆腔区、浅蓝色为塑性区)

Figure  3.  Characteristic dimensions of blasting zone formed by two kinds of charge(dark blue is the blasting zone and light blue is the plastic zone)

下载: 全尺寸图片 幻灯片

根据爆腔形状的模拟结果可以看出，柱形药包形成的爆腔表面较平滑，而等体积分成球形药包后形成的爆腔带有轻微的锯齿结构。该结果主要是因为，球形药包之间有弧形间隙，导致初始爆炸冲击波方向不同，并相互叠加不均匀，作用在土介质上形成锯齿结构。从爆腔尺寸上看，两种划分方式形成的爆腔形状基本一致，均近似为椭圆形，6 kg柱形药包等体积划分为8个球时，形成的爆腔和塑性区特征尺寸与柱形药包形成的爆腔和塑性区特征尺寸接近，离心率相差较小，离心率分别为0.380 0、0.708 9，所以按照图1(a)进行划分较合理，此时式(14)～(15)中，N=n。球形药包爆腔的长半轴和短半轴分别约为球形药包半径的10.4倍和5.1倍，这在实验得到的空腔半径约为装药半径的11倍^[11]结果范围之内；而塑性区长半轴约为爆腔长半轴的1.6倍，塑性区短半轴约为爆腔短半轴的2.85倍。

根据一维球形理论计算模型的叠加结果形成的爆腔包络线和一维柱形理论计算模型的爆腔半径进行对比如表1和图3所示。N=8时，一维柱形药包爆腔理论计算半径为215.7 mm，一维球形药包爆腔理论计算半径为189.5 mm。

由图3可以看出，8个球形药包近似替代形成的爆腔两侧为锯齿结构，数值模拟结果与图2(a)的结果基本一致，其爆腔短半轴比一维柱形爆腔理论半径略小，相差11.9%。这是因为一维柱形爆腔理论计算为无限长柱，轴向无冲击波能量传播，不考虑轴向能量损失，只考虑径向能量变化，所以对于有限长柱形药包存在轴向端头效应使能量向轴向传输，实际结果应该更接近球形叠加爆腔尺寸。

由计算结果可以看出，采用8个球形药包叠加的近似计算方式得到的结果与柱形药包模拟计算结果较一致，爆腔长半轴相差1%，短半轴相差5.25%，离心率相差4.29%。结果表明，可以用一维球形药包理论计算叠加的计算方式来计算二维柱形药包爆腔特征尺寸，能满足工程需要。

根据不同的计算方式得出爆腔特征尺寸计算结果见表2。

表  2  不同计算方式爆腔特征尺寸计算结果

Table  2.  Characteristic sizes of blasting cavity calculated by different calculation methods

计算方式	半径/mm	爆腔长半轴/mm	爆腔短半轴/mm	离心率
柱形模拟	38.84	530	200	0.377 4
8个球形替代模拟	47.89	500	190	0.380 0
一维柱形理论计算	38.84		215.7
8个球形叠加近似计算	47.89	524.7	189.5	0.361 2

下载: 导出CSV 

| 显示表格

3.   不同药量、不同长径比药包的计算结果分析

按图1(a)的柱形药包转球形药包划分方法，利用式(13)～(14)，分别对8、10、12 kg TNT长径比为10、12、15、17、20的情况进行数值模拟和球形药包叠加形成的爆腔、塑性区特征尺寸及误差对比分析，结果如图4～5所示。

图  4  不同药量、不同长径比情况下，爆腔特征尺寸理论计算结果与数值模拟结果的对比和相对误差分析

Figure  4.  Comparison between theoretical calculation results and numerical simulation ones for characteristic dimensions of blasting chambers with different charge mass and different length-to-diameter ratios, and relative error analysis

下载: 全尺寸图片 幻灯片

图  5  不同药量、不同长径比情况下，塑性区特征尺寸近似计算结果与数值模拟结果的对比以及相对误差分析

Figure  5.  Comparisons between theoretical calculation results and numerical simulation ones for characteristic sizes of plastic zones at different charge mass and different length-to-diameter ratios, and relative error analysis (left)

下载: 全尺寸图片 幻灯片

从图4～5可以看出，进行球形理论模型叠加后形成的爆腔与柱形药包数值模拟结果较接近，且理论计算结果小于模拟结果，误差趋势随长径比的增大而减小，与模拟结果比较：

（1）爆腔长半轴误差范围为1.30%～8.75%，短半轴误差范围为3.9%～12.2%；

（2）塑性区长半轴误差范围为1.7%～6.4%，短半轴误差范围为2.6%～7.3%。可以看出，两者长半轴误差较小，短半轴误差略大。

产生误差的原因主要有两方面：

（1）理论计算结果的爆腔尺寸未考虑其他药包径向和轴向的叠加效果，只考虑端头药包的爆腔及塑性区尺寸，且产生塑性区所需冲击波压力小于形成爆腔的冲击波压力，此时忽略药包径向叠加效果会使塑性区短半轴尺寸产生误差；

（2）炸药在土介质中爆炸后产生地震波，在自由界面进行反射，反射波对介质的拉伸会进一步扩大爆腔尺寸，而理论计算结果未考虑反射波的影响，结果较小。

同时，相同药量长径比越大，按照式(13)～(14)的球形药包叠加近似计算方法形成的爆腔特征尺寸(爆腔长半轴和短半轴)越接近柱形药包爆腔的模拟结果，当长径比较大的时候，球形药包划分个数变多，与土介质接触的侧面趋向于一条直线，土介质受到的径向应力分布的越均衡，更接近柱形药包爆炸的应力荷载分布。

4.   结　论

（1）一维球形装药爆腔预测准静态理论模型能够利用叠加的方法，近似计算柱形装药爆腔的特征尺寸。当球形药包数N=n、药柱长径比在10及以上时，该计算方法较合理，其爆腔尺寸和塑性区尺寸的计算结果均小于数值模拟结果，误差较小，在12.2%以内。因此，该方法可以用于爆腔的工程计算。

（2）随着柱形药包长径比的增大，球形叠加模型近似计算结果的精度增高。球形叠加模型近似计算方法适用于长径比较大的柱形药包爆腔尺寸及塑性区尺寸的计算，且划分方式简单易算，表达式中参数获取较方便，能够快速给出工程结果。