基于损伤-虚拟张拉裂纹模型的地下爆炸围岩破坏规律研究*

2022-06-14程树范曾亚武

爆炸与冲击 2022年5期

程树范，叶阳，曾亚武，高睿

（武汉大学土木建筑工程学院，湖北武汉 430072）

球形硐室三维力学性能优越，在军事工程中应用广泛。开挖形成的球形腔室，不需要复杂的支护即可作为隐蔽空间进行爆炸实验。作为爆炸危险物的贮藏点，硐室存在潜在的内爆炸风险，弄清硐室抗爆性能和爆炸后的围岩损伤规律，对硐室的防护及修复十分重要。

与填实爆炸相比，大型空腔的间隔作用有效规避了爆炸产物对腔壁的直接冲击，并削弱了爆炸冲击波的强度，有利于维持周边岩体的稳定性。李孝兰对空腔解耦效应进行了系统总结，认为随着空腔半径的增大，爆炸荷载最终能够实现解耦。楼沩涛等在硬岩场地中进行了球型空腔的间隔爆炸试验，结果表明，半径1.85 m 的球型硐室基本可以实现240 kg TNT 炸药的完全解耦；王占江则进行了微量装药的解耦模拟试验，也验证了空腔解耦的可行性。由于爆炸试验危险性高，且试验结果难以被准确地量测，因此数值模拟被广泛用于地下结构的抗爆设计，如：Wang 等采用数值模拟研究了多源爆炸条件下地下硐室的动态响应，认为应加强爆炸源中段的防护措施；熊益波等基于数值模拟技术对抗爆硐室可能存在的高风险区进行了预测，并制定了有针对性的加强方案。对于爆炸问题，在连续介质理论的基础上，针对岩体大变形和高应变率工况，研究人员有针对性地提出并改进了材料模型，使得模拟结果更加可靠，然而在模拟脆性开裂时仍存在明显的不足。而离散元方法虽然可以较好地模拟裂纹的形成与扩展，但并不适合大尺度的模拟。目前，部分研究基于“生死”单元，在有限元模型中通过单元失效模拟裂纹，取得了较好的模拟效果，但单元删除导致的质量不守恒和能量异常损失仍会降低数值模拟的可靠性，特别是当网格尺寸较大或破碎程度较严重时，该方法的模拟结果将严重失真。

为克服连续介质方法的不足，本文中在岩石HJC 模型的基础上，提出一种新型的损伤-虚拟裂纹模型，首先讨论该模型的适用性和尺寸效应，并基于该模型对硐室内爆炸引起的围岩损伤规律进行分析，最后通过与有限元方法的对比，讨论该模型及模拟方法的优越性。

1 数值模拟方法

1.1 岩石材料的HJC 模型

数值模拟过程中，材料模型的选择是否合理直接决定了数值模拟的可靠性，Holmquist 等在金属材料Johnson-Cook 模型基础上提出的HJC 模型综合考虑了静水压力、应变率和损伤对材料强度及力学性能的影响，是模拟脆性材料动态破坏的一种较理想的模型。HJC 模型有19 个独立的模型参数，分别为基本力学参数、、、ρ，极限面参数、、、，压力参数、、、、、、，应变率参数和损伤参数ε、、。HJC 模型虽然包含较多参数，但其物理意义明确，标定方法也已经较为成熟。

本文研究对象为我国西南地区完整性较好的硬质红砂岩，通过单轴及常规三轴压缩实验、巴西劈裂实验，得到自然状态下红砂岩的基础力学参数，见表1。

表1 红砂岩的基本力学参数Table 1 Basic mechanical parameters of red sandstone

表1 中单轴压缩实验对应的加载应变率为10s，为标定应变率参数，另外还进行了3 组不同应变率的单轴压缩实验和1 组分离式霍普金森压杆(split Hopkinson pressure bar，SHPB)实验，得到应变率为5×10、2×10、10和87.3 s时硬质砂岩的动态强度分别为79.18、83.57、85.13 和100.65 MPa。最后根据文献[18-19]中所建议的方法，得到HJC 模型参数，见表2。

表2 红砂岩HJC 模型参数Table 2 Parameters of the HJC model of red sandstone

采用表2 中的材料参数，基于有限元方法（FEM）在LS-DYNA 软件平台上进行了单轴压缩的数值模拟，以验证所标定参数的可靠性。考虑到数值计算的效率，根据对称性，将圆柱体试件简化为图1(a)所示的二维模型进行分析，模型尺寸为50 mm×100 mm，采用非结构化的Delaunry 网格划分，网格控制尺寸为2 mm，压力板与试件间采用基于罚函数的双向接触，接触摩擦因数为0.1。模拟得到的单轴压缩应力-应变曲线和破坏后的形态如图1(b)所示。

图1 单轴压缩数值模拟结果Fig. 1 Numerical simulation results of uniaxial compression

由图1(b)可知，单轴压缩数值模拟得到的应力-应变曲线与实验值一致性较好。虽然FEM 模型的最终破坏形态与实际状态存在一些差异，但均以剪切破坏为主，且模拟得到的破坏倾角与主裂纹倾角十分接近，对应的破坏机理是一致的。可见经过参数标定的HJC 模型在模拟岩石准静态单轴压缩时的结果是可靠的。

相应地，为验证其在动态破坏模拟时的有效性，还进行了SHPB 实验的数值模拟，建立的实验装置几何模型如图2 所示。

图2 中的入射杆和透射杆均为低合金钢材质，弹性模量为210 GPa，泊松比为0.27。入射杆、透射杆的长度分别为2.5、1.5 m，直径为50 mm；岩石试件高度和直径均为50 mm，杆件和试件的网格尺寸分别为2、5 mm。对于入射应力波的加载，采取时程荷载的方式来实现，即将实验测得的入射波直接施加于入射杆，因此数值模型中并不包括子弹。实验与数值模拟得到的透射波、入（反）射波以及由三波法得到的应力-应变曲线如图3 所示。

图2 SHPB 冲击压缩实验的数值模型Fig. 2 Numerical model of the SHPB impact test

图3 SHPB 冲击压缩数值实验结果Fig. 3 Numerical simulation results of the SHPB impact test

由于采用了时程荷载的方式施加入射波，图3(a)中入射波段完全重合，而模拟得到的反射波和透射波也具有较高的吻合度，数值模拟得到的加载平均应变率为90.36 s，动态抗压强度为103.22 MPa，相应的实验值分别为87.3 s和100.65 MPa，误差均小于5%，说明经标定的HJC 模型在模拟动态破坏时也具有较为理想的效果。

然而在本构关系方面，相较于压缩破坏时的分段考虑，HJC 模型在模拟岩石受拉破坏时采用的理想弹塑性本构并不适合描述岩石的脆性断裂。采用表2 中的参数，基于FEM 模拟了巴西劈裂实验，模拟方法与单轴压缩试验相同。得到的破坏形态及应力-应变曲线如图4 所示，与室内实验有较大的差异，试件最终表现为与单轴压缩相似的斜裂纹破坏。这是由于HJC 模型仅在受压时考虑了损伤对材料力学性能的影响，而损伤变量的定义又包含了张拉产生的塑性应变，在静水压力为拉的情况下，张拉产生的损伤会导致材料抗压强度快速下降，产生图4 所示的异常破坏形式，这是HJC 模型的一个重大不足。

图4 采用有限元方法进行的巴西劈裂数值模拟Fig. 4 Numerical test of Brazilian splitting test based on FEM

1.2 考虑张拉失效的虚拟裂纹模型

为克服HJC 模型在模拟张拉破坏时的不足，考虑引入虚拟裂纹来模拟张拉破坏，而不考虑材料本身的张拉屈服，即采用如图5(a)所示的节理单元连接代替FEM 中所采用的实体单元共节点连接。建模过程中，首先将单元离散化，此时节点在空间上重合，但各自独立，其后在相邻实体单元间插入一无厚度节理单元，通过节理单元与实体单元间的共节点连接，最后将离散单元重新组合为整体。节理单元采用张拉失效的双线性内聚力模型，内聚力单元受力后可以在法向产生独立的位移，而在切向不会发生变形，其法向应力-位移曲线如图5(b)所示，包括无损弹性和损伤软化两个阶段。

图5 张拉失效的双线性内聚力模型Fig. 5 Cohesive model with a tension failure rule

图5 中，σ为法向应力，为法向刚度，δ为法向相对位移。δ为法向的损伤起始位移，当δ＜δ时，内聚力单元完全弹性，当δ＞δ时，损伤演化开始，δ=δ为极限状态；σ=δ为法向最大应力；δ为法向的失效位移，当δ＞δ时，虚拟裂纹被激活。

根据图5(b)，损伤起始后，内聚力单元的本构关系可以表示为：

式中：为法向的损伤变量，可定义为：

由于节理单元在切向不会出现位移，剪切变形需要通过实体单元来完成，因此，在压缩状态下，本文模型将退化为FEM 模型。插入节理单元后，将引入3 个细观材料参数、δ和δ，其取值无法通过实验直接得到，需要通过巴西劈裂的数值模拟反演得到，即要求模拟得到的荷载-位移曲线以及最终的破坏模式均与实验一致。最终标定的细观参数为：=4.2 GPa/mm，δ=0.012 mm，δ=0.025 mm，对应的荷载-位移曲线如图6(a)所示，无论是应力-应变曲线，还是最终的破坏状态，均与实验高度吻合。根据数值模拟结果，内聚力单元对单轴压缩过程的模拟结果影响较小，插入内聚力单元后的应力-应变曲线和破坏形态如图6(b)所示。

图6 巴西劈裂及单轴压缩数值模拟Fig. 6 Numerical simulation of the Brazilian splitting and uniaxial compression based on the proposed method

1.3 模拟方法的可靠性

虚拟裂纹的引入，虽然解决了HJC 模型在模拟低静水压力时的不足，但虚拟裂纹的插入不可避免地会导致数值模型出现网格依赖性和尺寸效应。且插入节理单元后，数值模型的计算量将显著增加，对于大型工程尺度的模拟计算，难以实现毫米级的网格划分，增大网格尺寸后，标定的模型参数是否适用，还需要进行专门的讨论。

HJC 模型的材料参数本身与单元尺寸无关，在模型被放大后，理论上不需要重新标定，但内聚力单元的参数与网格尺寸密切相关，对于张拉破坏，由于破坏面较为稳定，可以引入一个比例因子来考虑尺寸效应，即在保持破坏应变不变的情况下，调整模量和位移来控制尺寸效应。当网格尺寸扩大倍后，对应的法向刚度参数保持不变，法向位移参数δ和δ则同步扩大倍，这样处理即可保持变形模量和破坏应变不变。为验证这一方法的可靠性，在保持网格形式一致的前提下，分别建立50 mm×100 mm、50 cm×100 cm、50 dm×100 dm 以及50 m×100 m 等4 个尺度的单轴压缩模型和单轴拉伸模型，对应的比例因子分别为1、10、100 和1 000。将单轴加载速率依次设置为2.5 mm/s、2.5 cm/s、2.5 dm/s 和2.5 m/s，以保证加载应变率一致，数值模拟得到的应力-应变曲线如图7 所示。

图7 不同尺度的单轴压缩和单轴拉伸应力-应变曲线Fig. 7 Stress-strain curves of uniaxial compression and uniaxial tension at different scales

由于本文模型的剪切破坏受节理单元的影响较小，在保持网格划分形式不变的基础上，随着网格尺寸的增大，材料强度和变形模量仅出现了小幅度的下降。在图7(a)中，当网格尺寸为2 cm 和2 dm 时，模型的强度为2 mm 网格时的94.55%和91.75%，减小幅度未超过10%，因此采用厘米和分米级的网格进行工程模拟也具有较高的可行性，且10%以内的误差可以作为安全储备，模拟结果偏于安全。当网格尺寸达到米的量级时，误差仍能保持在20%以内，模拟结果尚可用于定性分析和粗略估算。根据图7(b)，当网格尺度由毫米量级增加至米量级后，相应的抗拉强度增加仅为4.2%，说明本文模型通过比例因子法较好地解决了虚拟裂纹带来的网格尺寸效应，可适用于多种尺度网格。由于本文方法同时考虑了岩石材料在高静水压力状态下的应变率效应、塑性损伤和低静水压力状态下的脆性断裂，因此适合于大尺度的地下爆炸的模拟。

2 地下填实及空腔爆炸的模拟

2.1 算例及建模方法

根据文献[6]，硬岩硐室的完全解耦（围岩对爆炸荷载的响应为完全弹性）时的比例半径约为0.35 m/kg，那么当TNT 炸药装药半径=0.5 m（装药量863 kg）时，解耦半径约为3.3 m。据此分别建立填实爆炸模型，和空腔半径=1、2、3.5 和5 m 的4 个空腔爆炸模型。根据试算，岩体部分模型尺寸取为20 m×20 m，空气域范围为6 m×6 m，即可实现爆炸荷载的有效耦合且虚拟裂纹及损伤均不会超出模型区域，以=2 m为例，建立的数值模型如图8 所示。

图8 空腔爆炸数值模型(R=2 m)Fig. 8 Numerical model of cavity explosion in LS-DYNA (R=2 m)

数值模拟的建立在显式动力分析软件LS-DYNA 中完成，其中岩体部分包括实体单元和节理单元，采用基于Delaunry 三角划分的Lagrange 网格，空气域则采用结构化的Euler 网格划分，空气域为正方形，并通过关键字Initial_Volume_Fraction_Geometry 对炸药材料空间位置进行定义，以提高网格质量。岩体与空气重叠部分通过罚函数实现流固耦合，由关键字Constrained_Lagrange_In_Solid 控制，Lagrange 网格和Euler 的控制尺寸均为0.2 m，且爆炸产物侵蚀岩体时，不允许网格间出现穿透，岩体外侧为无反射边界。

LS-DYNA 软件自身拥有丰富的材料库，其中High_Explosive_Burn 材料被广泛用于凝聚相炸药的模拟，以TNT 炸药为例，其爆速=6 930 m/s，爆压=21 GPa。爆炸产物的扩容过程则通过附加JWL（Jones-Wilkins-Lee）状态方程来实现，其表达式为：

式中：、为线性爆炸系数，具有压强的量纲；ω、和为无量纲爆炸系数；为相对体积，即压强为时气体体积与初始体积之比；为炸药单位体积爆轰能量。

选取的TNT 炸药参数见表3。

表3 TNT 炸药爆轰产物JWL 参数Table 3 JWL EOS parameters of the TNT detonation product

炸药为中心起爆，爆炸过程空气可视为理想气体，其状态方程为：

式中：为空气压力；ρ/ρ为相对密度，其中ρ 为压缩后空气的密度，ρ为空气初始密度，取为1.29 kg/m；γ 为理想气体比热容，取为1.4；为单位体积内能，取为250 kJ/m。

在LS-DYNA 中，理想气体可用附加多项式状态方程的Null 材料模拟，标准状态方程为：

式中：=(ρ -ρ/ρ；～为多项式系数，当==0、γ -1、0 时，式(5)与式(4)完全等效。

2.2 爆炸荷载

对于地下填实爆炸，爆炸荷载以爆炸产物压力的形式直接作用于岩体，其峰值压力可达10 GPa 数量级；而对于空腔爆炸，由于空气的间隔作用，爆炸荷载多以空气冲击波的形式作用于岩体，数量级相对较低。填实和空腔爆炸在荷载形式上有较大区别，荷载的计算方法也有所差异。

对于填实爆炸，爆炸产物压力可以采用等熵膨胀法进行计算，亨利奇给出的耦合装药最大压力的计算方法为：

式中：为作用于岩体的最大压力；ρ为炸药密度，取为1 650 kg/m（TNT）；为炸药爆轰速度，取为6 930 m/s；γ 为等熵指数，取为3。

而对于空气冲击波压力，则可以采用最大扩散速度法进行计算，其表达式为：

式中：ρ为空气密度；为岩石内冲击波速，当冲击波强度较小时，可取为剪切波速1 264 m/s；为空气绝热指数，对于双原子分子可取为1.41；为冲击波压力增大系数，取值为0～20。

以=0.5， 5 m 的模型为例，由式(6) 计算得到的填实爆炸对应的孔壁最大压力为9.91 GPa；由式(7)计算得到的空腔爆炸最大压力范围为0～34.2 MPa。而本文数值模拟监测到的孔壁最大压力分别为8.45 GPa 和10.3 MPa，均在理论值附近，说明ALE 算法能够较为准确地模拟爆炸荷载作用。

2.3 围岩分区破坏规律

数值模拟结果见图9，地下填实爆炸过程中围岩的损伤在时间和空间上都具有很强的分区特征。

在时间上可以分为3 个阶段：第1 阶段为孔壁的外扩，此阶段内爆炸气体直接冲击岩体，产生极大的压力，远超岩石的屈服强度，岩石表现出一种类似于流体的力学性质，岩体向外移动形成爆炸空腔，并激发压缩波；第2 阶段为破碎区的形成，此阶段内压缩波强度大于岩石的屈服强度，岩体出现大的塑性损伤，岩体完全破碎；第3 阶段为裂隙区的形成和扩展，此阶段内压缩波压力小于岩体的抗压强度，但受到泊松效应影响，环向的拉应力仍可诱发张拉裂纹的扩展。随着压缩波压力下降至不足以诱发裂纹继续扩展后，爆炸过程结束。相应地，根据破坏机制，在空间上则会形成扩容空腔（见图9(a)）、压缩破坏区和张拉破坏区（见图9(c)）。

根据第1.2 节内容，内聚力单元的失效可以作为裂纹产生和扩展的可靠判据，而实体单元的损伤则可以作为岩石屈服的判据。相应地，对于地下爆炸，两者可用于裂隙区和破碎区范围的确定。考虑到装药量对爆炸效果的影响，根据霍普金斯比例定律，引入比例半径来描述地下爆炸的分区破坏规律。比例半径的定义为：

式中：为爆心距；为装药质量，对于本文模型，可取=863 kg。

如图9(c)所示，耦合装药时破碎区平均半径为2.0 m，是炮孔半径的4 倍；考虑装药尺寸后的比例半径为0.26 m/kg；而裂隙区半径为5.0 m，是炮孔半径的10 倍，对应的比例半径为0.65 m/kg。

图9 地下填实爆炸的分区破坏规律Fig. 9 Zonal failure law of the underground coupled explosion

2.4 硐室内爆炸的空腔解耦

数值模拟得到的空腔爆炸模型，其最终的破坏形式如图10 所示。

对于空腔爆炸，破碎区的厚度随着空腔尺寸的增大而减小，说明空腔效应可以减小作用于孔壁上的爆炸压力，从而提高硐室的抗爆性能。当硐室半径达到预设的解耦尺寸（=3.5 m）时，围岩表面仅出现零星的实体单元损伤，而壁面附近却分布有密集的短裂纹。分析认为，当爆炸冲击波在硐室表面反射时会产生反射，反射波与后续的冲击波叠加使爆炸压力被加强，并转化为张拉应力波，由于岩体的抗拉强度远小于抗压强度，张拉作用仍能引起表层的局部破坏，这一现象在文献[10]中也被观测和验证，这部分由反射波叠加形成的裂纹虽然十分密集但没有扩展的趋势，因此仍将其归于裂隙区。当硐室半径达到5 m 后，裂隙区完全消失，硐室处于完全弹性状态，据此认为：对于硬质砂岩硐室，比例半径为0.52 m/kg时，可以实现爆炸荷载的完全解耦。

3 讨论

对本文损伤-虚拟裂纹模型与传统FEM 模型的差异与优劣，进行简单讨论。

在图10(a)所示的耦合装药模型的基础上，通过调整装药半径，可以将装药形式转化为不耦合装药，从而减小围岩损伤，当装药半径调整为0.12 m 时，对应的损伤及虚拟裂纹如图11(a)～(b)所示，围岩中破碎区范围较小，且出现了多条细长裂纹，形成半径2 m 的裂隙区。

图10 空腔爆炸模型最终的破坏形态Fig. 10 Final failure characteristics of the cavity explosion models

图11(c)为相同条件（除无黏聚力单元外，其他条件均与本文模型相同）的有限元方法模拟的损伤云图。对比图11(b)～(c)可知：FEM 模型计算得到的损伤具有很强的对称性，损伤云图近似呈环状，而本文的损伤-虚拟裂纹模型计算得到的损伤演化受裂纹扩展的影响较明显，裂纹延伸处的损伤更加严重，对应云图呈现放射状，体现了裂纹对损伤演化过程的导向作用，虚拟裂纹的随机性即可反映出岩石的各向异性特征。由于节理单元失效消耗了一部分能量，有限元总损伤面积和损伤程度相对较小，但由于裂纹分布不均匀，实际的损伤演化范围较有限元模型略大。综上所述，损伤-虚拟裂纹模型能够考虑到岩体裂纹对损伤分布的导向作用，综合考虑损伤及裂纹效应后的模拟结果在理论上应该更加真实。