岩石爆破损伤模型的比选与改进

胡英国 ，卢文波 ，陈 明 ，严 鹏 ，周创兵

（1. 武汉大学 水资源与水电工程科学国家重点实验室，武汉 430072；2. 武汉大学 水工岩石力学教育部重点实验室，武汉 430072）

1 引 言

岩体爆破是大型水利、采矿工程开挖必不可少的施工手段，爆破开挖将对围岩造成一定程度的损伤，影响工程安全，因此，研究爆破损伤区特性，从而有效地控制围岩的损伤范围，具有重要的工程意义。岩体的爆破损伤是一个复杂的动态演化过程，为了研究探索这一过程，研究者们相继建立了相关的岩体爆破损伤模型。

美国Sandia国家实验室早在1966年就开始了岩石爆破损伤模型的研究工作[1]，主要方法是将岩石的动态断裂作为一个连续的损伤累积过程来处理，其基本点是建立损伤变量与岩石内微裂纹密度的关系，并预测在爆炸载荷作用下岩石的损伤和破坏过程。Grady和Kipp[2]提出了岩石爆破各向同性损伤模型，即GK模型，该模型采用一个标量描述被拉应力激活的钱币状裂纹所引起的岩石刚度的劣化，同时假定这些裂纹数服从双参数的 Weibull分布，他们采用该模型模拟爆炸载荷作用下油页岩的动态断裂和破碎，并根据能量平衡准则得到了与应变率有关的碎块平均尺寸表达式。Taylor等[3]引进O 'connell、Budianshy[4]的关于有效体积模量和泊松比与裂纹密度的关系表达式以及 Grady 给出的碎块尺寸表达式，建立了损伤变量与裂纹密度之间的关系式，并将损伤变量以率形式耦合到动态本构方程中，该TCK模型可以预报岩石在体积拉伸载荷下的动态响应。Kuszmaul[5]在以上两模型的基础上提出了KUS模型，该模型考虑了高密度微裂纹下的荫屏效应，即微裂纹周围产生应力释放的材料能够重叠，在裂纹的激活率中考虑了损伤引起的减少。Thorne等[6]在前人的基础上考虑了激活裂纹数可能引起岩石体积的变化，并通过采用不同的损伤变量定义，考虑了模型在大裂纹密度条件下的适应性，建立了 Thorne模型。Yang、Liu等[7－8]对以上模型在裂纹密度的分布及损伤变量的定义方面进行了修正，认为只有在体积应变大于某一临界体积应变后裂纹才能扩展，并考虑作用时间对裂纹密度的影响，在定义损伤变量时引入了断裂概率的概念。随着计算机硬件技术的发展，涌现出新的研究岩体爆破损伤过程的方法，唐春安等[9]不引进裂纹密度，采用累计拉应变和极限应变的比例关系来反映岩体的损伤程度，开发了用于研究岩体断裂破坏过程的RFPA软件。朱哲明等[10]将应力判据嵌入AUTODYN，并对应力更新算法进行修正，研究并揭示了岩体爆破损伤区形成的力学机制。王志亮等[11]基于 LSDYNA自定义接口将TCK模型导入LS-DYNA并计算，详细研究了爆破拉伸损伤区的特性；马国伟等[12]采用同样的方法将Johnson- Holmquist模型嵌入 LS-DYNA，研究了爆破裂纹扩展的影响因素，并对爆破损伤控制提出建议。

采用损伤力学来研究问题时，其主要步骤为：首先定义一个合适的损伤变量，然后根据外荷载确定研究对象在荷载作用下的损伤演化方程和考虑损伤的本构关系。因此，反映岩体损伤程度的损伤变量的定义是基础，正确、合理的损伤变量不仅能够使要研究的问题简单明了，而且计算结果也更加准确。已有的爆破损伤模型的损伤变量定义方式多样，表1给出了5种应用较为广泛的爆破损伤模型对应的损伤变量的表达方式及相关内变量定义的对比（严格意义上来说RFPA不是一种爆破损伤模型，但其中损伤变量的定义方法值得关注和借鉴，在这里为统一称之为RFPA模型）。

从表1中可以看出，TCK模型与KUS模型的损伤变量表达式完全相同，只是在裂纹密度的表达中考虑了损伤的影响，而 THRONE模型的裂纹密度的表达式与TCK、KUS相似，但损伤变量的表达式差别明显，YANG等的损伤模型的裂纹密度和损伤变量表达式与前三者均有差别，RFPA模型的损伤变量的定义方法与其他几种模型区别最为明显，没有建立裂纹密度与损伤变量的关系，其表达式相对简单。上述爆破损伤模型虽然在相关的研究中取得了较好的效果，但其适用性和精确性并没有基于具体的工程实例进行系统地比较和验证，本次研究中，将基于具体的工程实例，对上述5种爆破损伤模型进行对比计算和验证。

表1 不同损伤模型相关变量的定义对比Table 1 Definitions of several variables for different damage models

2 岩体爆破损伤模型的比选

由弹性力学的基本理论可知，声波在岩体中的传播速度与岩体的物理力学性质密切相关，下面将给出通过声波测试判据得到岩体损伤区的理论依据，岩体中纵波传播速度可表示为

式中：VP为岩体声波速度；ρ为岩体密度；K为体积模量；ν为材料的泊松比。定义岩体的声波波速降低率η为：

采用η表示有效体积模量，若近似认为密度和泊松比不变，根据式（2）可得式（3）：

对比式（3）、（4）可以看出，声波波速与爆破损伤模型描述爆破开挖的损伤效应是等价的，只是前者从宏观声波速度衰减来体现，而后者采用与微观的裂纹密度密切相关的损伤变量来表征。

根据损伤变量式（3）、（4）可得爆破损伤模型的损伤变量D与声波降低率η的对应关系为

根据《水工建筑物岩石基础开挖工程技术规范》[13]规定，当波速降低率大于10%时，即判定岩体受到爆破损伤破坏，代入式（5）计算可得爆破损伤模型的损伤变量D的损伤阈值为0.19。

为方便进行爆破损伤模型的比选计算，宜选择应力状态较为简单的工程实例以减少其他因素的干扰，本次研究选择夏祥等[14]、李海波等[15]对广东岭澳核电站基岩爆破损伤区进行声波检测的实测结果为工程依据，文献[14－15]提供了广东岭澳核电站的爆前、爆后的声波监测资料，并给出了围岩体损伤范围。图1、2为声波检测测点布置示意图以及围岩体损伤区示意图。

图1 声波检测布置示意图Fig.1 Arrangement of acoustic testing

图2 损伤区实测示意图Fig.2 Measured depth of damage zone in rock

根据相关的文献资料，建立如图3所示的模型，模型中圆柱直径为 10 m，高为 8 m。装药半径为25 mm，同时计算过程中岩体参数以及爆破损伤模型的相关参数如表2、3所示[3,14]。计算过程中，通过在孔壁上施加荷载曲线以模拟爆炸荷载的作用，爆炸荷载随时间变化历程参考文献[16]的方法确定，其曲线形状如图4所示，边界条件设为无反射边界以减小应力波反射的影响，反映爆源周围岩体的实际受力状态。

图3 模型示意图Fig.3 Model sketch

表2 模型相关物理参数Table 2 Pysical parameters of model

表3 计算中损伤模型对应的损伤参数Table 3 Damage parameters of damage models

图4 爆炸荷载随时间变化曲线[16]Fig.4 Curves of blasting load versus time[16]

要实现以上爆破损伤模型的对比计算，需要借助 LS-DYNA的自定义接口，将其导入 LS-DYNA进行计算。实现LS-DYNA的二次开发的关键是用FORTRAN 语言编写程序描述该模型建立的数学过程的子程序，在每个积分步，LS-DYNA主程序与该子程序进行数据传递。LS-DYNA中自定义方法在其用户手册[17－18]中有详细介绍，这里不再赘述。

基于以上方法，实现上述5种爆破损伤模型，采用自定义接口将其成功导入LS-DYNA进行对比计算，图5分别给出了5种爆破损伤模型对应的损伤区云图。

图5 不同损伤模型损伤区域云图Fig.5 Nephogram of damage zone of different models

计算结果表明，从损伤区的分布来看，几种爆破损伤模型的分布规律均与实测值大致吻合，围岩中形成柱状损伤区，损伤区的损伤程度随着爆心距的增加而降低，炮孔附近岩体的损伤变量接近 1，由于炸药爆炸瞬间产生的巨大的爆轰压力，形成压致粉碎区，岩体完全破坏，随着远离炮孔壁，冲击波衰减为应力波，不能直接压碎岩体，但岩体切向的拉应力超过其抗拉强度，围岩中出现拉裂纹或原有裂纹进一步扩张，出现中等损伤程度的裂隙区，随着爆炸荷载的进一步衰减，裂纹停止扩张，形成未损伤区。TCK模型、KUS模型和THRONE模型对应的损伤区外沿的整体形状非常类似，只是不同位置对应的损伤程度以及损伤范围不同；TCK模型的损伤范围最大，而 THRONE的损伤范围最小；YANG等建议的损伤模型损伤区的外沿与前3个模型偏差较大，由于RFPA模型损伤变量的定义方法与以上4种明显不同，从计算结果可以看出，RFPA模型对应的损伤区，无论是损伤区外沿上还是损伤区内岩体对应的损伤程度与其他4种损伤模型区别相对明显，由此可见，不同的损伤变量定义方法计算得到的损伤区的空间分布特征存在明显区别。

为进一步比较不同损伤模型模拟爆破损伤区的准确性，表4给出了不同模型算例中损伤区的最大水平半径与孔底损伤深度。

表4 不同损伤模型的损伤范围对比Table 4 Damage area comparison between different models

表4数据表明，对于损伤最大水平半径，最接近实测值的是KUS模型和RFPA模型，其中TCK模型、KUS模型、YANG模型偏大，最大误差为0.4 m，对于孔底的损伤深度，除了PFRA外，其余几种损伤模型的孔底损伤深度均偏小，其中误差最大的是YANG等的损伤模型，其误差达到了0.58 m，最接近实测资料的依然是KUS模型和RFPA模型。因此，计算结果表明，5种爆破损伤模型中，KUS模型和RFPA模型对应的损伤区范围最接近实测结果。

3 拉压损伤模型的建立与验证

以上爆破损伤模型在考虑损伤效应时，均将拉应变与损伤变量建立直接的关系，这些损伤模型可以较好地反映岩体爆破产生的拉伸损伤，但不能有效地体现压损伤的特点。在炸药爆炸瞬间，巨大的爆轰压力在孔壁周围形成一定范围的压致粉碎区，这一点在上述模型中没有考虑。因此，有必要建立能同时反应拉伸、压缩损伤的损伤模型，使其具有更广泛的应用价值。根据上文的计算结果，KUS模型计算的损伤范围的精确性相对较高，本次计算将在对 KUS模型进行改进的基础上，考虑引入压损伤，建立拉压损伤模型，下面将给出拉压损伤模型的详细数学推导过程。

根据 Taylor[3]等对损伤体积模量与裂纹密度的相关关系的计算表达式：

Grady等[2]认为，裂纹密度是裂纹影响区岩石体积与岩石总体积之比，激活的裂纹数服从体积拉伸应变的双参数Weibull分布，即

式中：γ为随机分布参数；N为单位体积裂纹数；a为微裂纹平均半径。

爆炸应力波作用下的微裂纹平均半径及裂纹密度可由式（8）、（9）决定：

Kuszmaul[5]考虑了高密度裂纹周围应力释放区的材料重叠，在裂纹的激活率中考虑了损伤引起的减少，将裂纹密度表达式修正为

式中：Dt为拉伸损伤变量。

考虑损伤效应的有效泊松比的表达式为

式中：θ为材料参数。

在体积拉伸状态下，拉伸损伤变量 Dt可由式（12）表达：

在压缩、剪切状态下，岩体服从各向同性的弹塑性本构关系，其屈服条件为

对于压缩损伤变量的考虑，引入 RDA模型对岩体中压缩损伤的处理方法，基于 RDA模型的应变率效应耦合原则，可得出在体积压缩下压缩损伤变量的表达式[18]为

由于模型中同时包含拉伸损伤变量 Dt和压缩损伤变量Dc，在表征岩体的损伤程度时，需要进行统一，在此定义损伤变量D，考虑损伤的最不利因素，D由下式决定：

得到损伤变量 D后，已有的爆破损伤模型（TCK、KUS等）的处理方法将所有宏观弹性常数均按一定比例减小，如：

以上处理方法存在一定的缺陷，它忽略了岩体的弹性常数之间需要满足的固定关系，通常情况下，岩体的宏观弹性常数应满足式（19）、（20）：

本次研究中，对宏观弹性常数的确定方法作如下修正，损伤体积模量和有效泊松比沿用已有的爆破损伤模型（TCK、KUS等）的处理方法，其他弹性参数则根据上式中的固定关系获得，改进后模型中的材料参数关系式如下所示：

其余宏观弹性常数由式（23）、（24）确定：

最后，爆破荷载下岩石的记录损伤效应的本构关系表达式可由以下增量型的虎克定律表示为

式中：σij为应力张量；εkk为体积应变张量；eij为偏应变张量；δij为Kronecker 符号。

通过以上推导，在对KUS模型进行修正的基础上，考虑压缩损伤建立了能够同时反映拉伸、压缩损伤的拉压损伤破坏模型。将以上过程通过FORTRAN 编程实现，并成功导入LS-DYNA自定义接口，这样就可以采用该模型进行计算并验证其有效性。将以上改进的爆破拉压损伤模型采用FORTRAN子程序实现，并成功导入LS-DYNA，对广东岭澳核电站基岩爆破损伤效应进行计算，得到相关的总体损伤范围与局部损伤范围如图6所示。

图6(a)表明，在KUS模型的基础上进行修正的拉压损伤模型最后形成的损伤区最大水平半径为6.87 m，孔底损伤深度为2.10 m，相比原模型，无论是损伤区的分布还是范围都与实测值更为接近。图 6(b)给出了损伤区压损图，压损区的范围约为0.45 m，在总体损伤区中所占比例很小，岩体在爆破过程中，孔壁岩体受到强烈的冲击波作用，超过岩体的抗压强度，形成压致损伤区，但由于爆炸荷载衰减很快，这个过程很短，因此，形成的压碎区范围很小。

图6 拉压损伤模型损伤区Fig.6 Damage zones of tension-compression damage model

该模型由KUS模型的基础改进而来，因此，图7给出了损伤区实测值、拉压损伤模型以及KUS模型三者对应的损伤区对比示意图，以便进一步验证拉压爆破损伤模型的有效性。

图7 损伤区最大水平半径对比Fig.7 Comparison of maximum horizontal radius of damage zone

相比KUS模型，拉压损伤模型考虑了近区的压缩损伤效应，计算结果表明，拉压损伤模型与实测值的吻合性明显要好，无论是最大损伤的水平半径还是孔底的损伤深度均更加接近实测值，这表明本次研究中建立的拉压损伤模型能够更好地描述岩体中的爆破损伤效应。

4 结 论

（1）基于岭澳核电站爆破损伤区检测的实测资料，对5种爆破损伤模型进行比选，数值计算结果表明，KUS模型和 RFPA 模型与实测结果较为接近。

（2）根据比选结果，选取KUS模型，引进RDA的压缩损伤的考虑方法，修正宏观弹性常数的确定方法，建立新的拉压损伤模型，并成功地导入LS-DYNA，进行计算和验证。

（3）基于具体的工程实例对拉压损伤模型的精确性进行验证，结果表明，改进的拉压损伤模型可更有效地描述爆破开挖扰动下岩体的爆破损伤效应。

但在实际工程中， 爆破近区岩体的复杂本构关系难以精确模拟，爆破近区岩体的复杂的破坏类型的精确描述尚需要进一步探索研究，本次计算基于各向同性的假定，没有考虑实际工程中岩体存在的软弱结构面的影响，对于实际工程岩体的各向异性特征没有体现。总体而言，本次计算中采用的拉压损伤模型对工程研究具有一定的参考价值和进步意义，但还需进一步修正和完善。

[1]杨小林, 王树仁. 岩石爆破损伤模型及评述[J]. 工程爆破, 1999, 5(3): 71－75.YANG Xiao-lin, WANG Shu-ren. Rock blast damage models and comments[J]. Engineering Blasting, 1999,5(3): 71－75.

[2]GRADY D E, KIPP M E. Continuum modelling of explosive fracture in oil shale[J]. International Journal of Rock Mechanics and Mining Sciences &Geomechanics Abstracts, 1980, 17(3): 147－157.

[3]TAYLOR L M, CHEN E P, KUSZMAUL J S.Microcrack-induced damage accumulation in brittle rock under dynamic loading[J]. Computer Methods in Applied Mechanics and Engineering, 1986, 55(3): 301－320.

[4]BUDIANSKY B, O′CONNELL R J. Elastic moduli of a cracked solid[J]. International Journal of Solids and Structures, 1976, 12(2): 81－97.

[5]KUSZMAUL J S. A new constitution model for fragmentation of rock under dynamic loading[C]//Proceedings of the 2ndInternational Symposium on Rock Fragmentation by Blasting. Canada: Keystone, 1987: 412－423.

[6]THORNE B J, HOMMERT P J, BROWN B.Experimental and computational investigation of the fundamental mechanisms of cratering[C]//Proceedings of the 3rd International Symposium on Rock Fragmentation by Blasting. Brisbane, Australia: [s. n.], 1990: 117－122.

[7]YANG R, BRWDEN W F, KATSABANIS P D. A new constitutive model for blast damage[J]. International Journal of Rock Mechanics and Mining Sciences &Geomechanics Abstracts, 1996, 33(3): 245－254.

[8]LIU Li-qing, KATSABANIS P D. Development of a continuum damage model for blasting analysis[J].International Journal of Rock Mechanics and Mining Sciences, 1997, 34(2): 217－231.

[9]TANG C A. Numerical simulation of progressive rock failure and associated seismicity[J]. International Journal of Rock Mechanics and Mining Sciences, 1997,34(2): 249－261.

[10]ZHU Zhe-ming, XIE He-ping, MOHANTY B. Numerical investigation of blasting-induced damage in cylindrical rocks[J]. International Journal of Rock Mechanics and Mining Sciences, 2008, 45(2): 111－121.

[11]WANG Zhi-liang, LI Yong-chi, WANG J G. Numerical simulation of tensile damage and blast crater in brittle rock due to underground explosion[J]. International Journal of Rock Mechanics and Mining Sciences, 2007,44(5): 730－738.

[12]MA G W, AN XM. Numerical simulation of blasting-induced rock fractures[J]. International Journal of Rock Mechanics and Mining Sciences, 2008, 45(6):966－975.

[13]中国水利学会施工专业委员会爆破学组. SL47－94 水工建筑物岩石基础开挖工程施工技术规范[S]. 北京:水利电力出版社, 1995.

[14]夏祥, 李俊如, 李海波, 等. 广东岭澳核电站爆破开挖岩体损伤特征研究[J]. 岩石力学与工程学报, 2007,26(12): 2510－2516.XIA Xiang, LI Jun-ru, LI Hai-bo, et al. Study of damage characteristics of rock mass under blasting load in Ling’ao nuclear power station, Guangdong province[J].Chinese Journal of Rock Mechanics and Engineering,2007, 26(12): 2510－2516.

[15]LI Hai-bo, XIA Xiang, LI Jian-chun. Rock damage control in bedrock blasting excavation for a nuclear power plant[J]. International Journal of Rock Mechanics and Mining Sciences, 2011, 48(2): 210－218.

[16]卢文波, 杨建华, 陈明, 等. 深埋隧洞岩体开挖瞬态卸荷机制及等效数值模拟[J]. 岩石力学与工程学报, 2011,30(6): 1090－1096.LU Wen-bo, YANG Jian-hua, CHEN Ming, et al.Mechanism and equivalent numerical simulation of transient release of excavation load for deep tunnel[J].Chinese Journal of Rock Mechanics and Engineering,2011, 30(6): 1090－1096.

[17]LSTC. LS-DYNA theoretical manual[M]. California:Livermore Software Technology Corporation, 2003.

[18]LSTC. LS-DYNA keyword user’s manual[M]. California:Livermore Software Technology Corporation, 2003.

[19]FURLONG J R, ALME M L, DAVIS J F. Modeling the dynamic load/unload behavior of ceramics under impact loading RDA-TR-00.0-0001[R]. Arlington: R&D Associates, 1990.