孔隙水压力及炸药埋深对堤坝爆炸效应的影响分析

张智超，陈育民，刘汉龙，王维国

（1.河海大学 岩土力学与堤坝工程教育部重点实验室，南京 210098；2.河海大学 土木与交通学院，南京 210098）

1 引 言

我国水利水电工程众多，许多重要的水工建筑位于人口密集的地区，三峡水利枢纽工程的拦河大坝高达185 m，最高蓄水水位可达175 m，南水北调工程有着近3000 km长的引水渠道，这些都不可避免地涉及到堤坝的建设和运营。加之2011年“中央一号文件”关于加快水利改革发展的决定，更是将堤坝的兴修和维护推向了新的高潮。然而，正是由于水利事业的发展与繁荣，其在给人民群众带来巨大的经济效益的同时，也往往成为战争或者恐怖分子的重点打击对象。一旦此类堤坝工程在爆炸等袭击下发生破坏乃至溃坝，将会给人民的生命财产造成巨大的损失。从另一方面来说，人为地利用爆炸荷载对堰塞坝进行破坏，往往还能起到事半功倍的效果，汶川地震中形成的多个堰塞坝，就是通过爆破拆除而成功地化险为夷，避免了坝体突然溃决而导致重大洪灾的发生，保证了下游百姓的生命和财产安全[1－2]。因此，开展爆炸荷载下堤坝破坏形态的研究，具有重要的实用意义。

岩土体在爆炸荷载下的破坏形态最直观地表现为爆炸弹坑。以往的爆炸弹坑研究大多集中于水平场地，而专门针对堤坝爆炸弹坑的研究尚不多见。Schmidt和 Holsapple[3]通过理论分析和试验验证认为，在离心机中模拟爆炸是采用室内小模型模拟原型的最为有效的方法，论证了模型与原型在爆炸条件下的相似关系，并就爆炸弹坑问题进行了一系列研究；穆朝民等[4－5]进行了自由场地爆炸弹坑的试验与数值分析；Sauseville[6]、Zimmie 等[7－8]对堤坝爆炸弹坑及爆炸后的稳定性等进行了离心机试验模拟；刘军等[9]基于 LS-DYNA软件模拟了土石坝在爆炸荷载作用下的动力响应，但由于其施加的爆炸荷载为TM5-585-1手册确定的经验时程曲线，对触地爆炸下空气冲击波引起的感生地冲击未作考虑，因此，与实际情况有所偏差，特别是在近地表区域，直接地冲击和感生地冲击将发生复杂的叠加和混合，对分析模拟结果的影响更不宜忽略[10]。

影响堤坝爆炸破坏效应的因素有很多[11]，本文首先选取炸药埋深进行分析；同时，筑坝材料本身的性质也不容忽略，例如，湿土较之干土在爆炸荷载下的动力特性就有很大区别，动载作用下含水土体的孔隙水压力发展会对土体结构的动力响应产生巨大影响，这种孔压效应在以往的爆炸弹坑研究中尚不多见。因此，本文在LS-DYNA软件框架内，建立空气、炸药与均质堤坝3种物质的模型，综合考虑了爆炸的直接地冲击与感生地冲击，利用多物质ALE法，分别就影响爆炸弹坑形成的两个重要因素——炸药埋深和土体孔隙水压力进行模拟分析，通过观察后处理中的物质流动变形，得到堤坝在不同爆炸工况下的弹坑破坏形态，以期能够对水利工程中的堤坝安全设计以及爆破拆除堰塞坝等工程应用提供参考。

2 算法和材料模型

2.1 多物质ALE算法

ALE方法是在材料域和空间域外引入参考域，在参考网格上进行控制方程的求解。它综合了拉格朗日方法与欧拉方法的优点，既解决了拉格朗日方法下材料的严重变形，又克服了欧拉方法下移动边界引起的复杂性问题，为爆炸冲击波与结构的相互作用问题提供了较好的解决方案。

ALE方法分为单物质ALE方法和多物质ALE方法。多物质ALE方法容许在一个网格中包含多种物质材料，通过跟踪每种材料的边界，在相应的单元中进行物质交换和输送[12]。爆炸荷载下的堤坝动力响应问题涉及到土体、空气和炸药多种物质的大变形问题，需要同时处理界面大变形和网格大变形，而多物质 ALE方法的以下几个特点就很好地解决了这些问题：

（1）允许一个网格内含有多种物质，物质界面可以穿过网格，即炸药爆炸产物、土体和空气能够相互侵入原来彼此所占据的空间，形成新的物质分布形态。因此，可以通过物质的流动状态来判断弹坑的形成过程及最终形态，而不需要通过网格变形或者单元删除的方法。

（2）用重分重映处理网格的大变形，可以避免爆炸作用下单元网格畸变可能导致的计算中断。

（3）用level set方法、VOF方法等界面追踪方法处理界面的大变形，克服了欧拉法难以精确描述移动边界的缺陷。

因此，利用多物质ALE法能够较为合理地对堤坝爆炸弹坑进行模拟分析。

2.2 土体材料模型

数值计算的核心是本构模型，其对计算结果的准确性有重要影响。LS-DYNA中的土体本构模型[13]MAT_FHWA_SOIL可以模拟诸多材料特性，如爆炸冲击加载下的应变率效应、应变硬（软）化、考虑孔隙效应的弹性本构、修正Mohr-Coulomb屈服面以及孔隙水压力效应等，是较为合适的计算土体爆炸响应的模型。它的优越性已在文献[14]中充分验证，因此，本文选取其来描述均质土坝筑坝材料在爆炸荷载下的力学行为。

为了确保数值计算的有效性和稳定性，当剪应力较小时，该模型将标准Mohr-Coulomb屈服面修正为一个光滑的曲面，并且垂直于压力轴，其表达式为

式中：P为压力；φ为内摩擦角；J2为偏应力张量第二不变量；K(θ)为张量平面角的函数；c为黏聚力；AHYP为决定修正后的 Mohr-Coulomb屈服面和标准的Mohr-Coulomb屈服面相似程度的参数。

当AHYP=0时，式（1）表示的是标准Mohr-Coulomb屈服面，当AHYP采用较大值时，修正后的屈服面明显偏离标准Mohr-Coulomb屈服面，对于数值模拟来讲，AHYP的取值应该小于ccotφ，一般按照式（2）进行选取：

图1为标准Mohr-Coulomb屈服面与修正后的屈服面对比。除了在低应力区外，二者几乎一致。

图1 标准与修正后的Mohr-Coulomb屈服面对比Fig.1 Standard and modified Mohr-Coulomb yield surfaces

同时，为了令屈服面在低围压状态下的形状呈现为更符合实际的三角形，将标准库伦函数K()θ改进为[15－16]

含水土体的孔压发展对计算结果具有重要影响。孔压的增长会降低有效应力，并由此导致抗剪强度的下降。对于非饱和土，当气体体积在加载过程中被压缩到0时，孔隙水压力将开始增加，使得土骨架的有效应力降低。模型采用以下关系式来计算孔隙水压力u：

式中：εv为体积应变；Ksk为孔压系数[14]；ncur为当前孔隙率；D2为气体孔隙坍塌前控制孔隙水压力的材料常数。

参数Ksk影响空气孔隙坍塌后的孔压-体积应变曲线曲率。当D2相对于Ksk非常大时，在体积应变达到空气孔隙体积之前，孔压几乎为0；但当D2降低时，孔压开始增加。D2通过Skempton法中计算孔隙水压力的参数B得到：

式中：K为体积模量；n为孔隙率；s为饱和度。对于饱和土体，D2=0[17]，于是式（4）简化为

因此，参数Ksk直接影响土体孔压的发展，进而影响到土体有效应力的变化。本文将通过调整参数Ksk的大小来描述孔压上升以及有效应力下降的程度，以此来研究孔压效应对爆炸弹坑的影响。筑坝材料参数值将在各个算例中具体给出。

2.3 炸药材料模型及参数

采用高能炸药燃烧模型和 JWL状态方程来模拟所使用的TNT炸药[13]。JWL状态方程能够精确描述爆炸过程中爆轰产物的压力、体积、能量特性，其关系式为

式中：A、B、R1、R2、ω为材料参数；E01为爆轰产物单位体积的内能；V为单位体积装药产生的爆轰产物的体积；P为爆炸产生的压力。炸药材料参数取值分别为[18]：密度为1.63 g/cm3，爆轰速度为6930 m/s，CJ压力为21 GPa，A=371 GPa，B=3.23 GPa，R1=4.15，R2=0.95，ω=0.38，E01=8 GPa，V=1。

2.4 空气材料模型及参数

采用空材料（NULL）模型和线性多项式状态方程表示空气的本构关系[13]。将空气介质简化为非黏性理想气体，假设冲击波的膨胀为绝热过程，其线性多项式状态方程可简化为

式中：ρ为空气密度；γ为绝热指数；E为单位初始体积的内能。参数取值分别为[13]ρ=1.25 kg/m3，γ=1.4，E=0.25 MPa。

3 数值模拟方案

基于上述的土体本构模型，利用多物质 ALE法，分别就不同炸药埋深和不同孔压上升程度的堤坝弹坑效应进行数值模拟，通过多个模拟结果的比较分析，探讨不同爆炸工况和不同筑坝材料特性下的堤坝爆炸成坑效应。

按照文献[6]中的试验方案建立堤坝模型，坝高为4 m，坝底宽为20 m，坝顶宽为4 m，左、右坡度比都为 1:2。为了减少计算量，将模型简化成厚度为一个单元的准二维模型，并根据对称性，只建立了半模型（见图2），在后处理中再映射成为整体模型。共划分为8668个节点，8578个单元。设计药包尺寸为0.4 m×0.4 m×0.1 m，质量为26 kg。药包位置在坝体中心线上变化。建模时，首先只划分空气和土体单元，然后再根据需要，在计算前将某坐标范围内原来的单元材料属性修改为炸药单元，以此方便地设置不同的炸药埋深。坝顶和坝坡与空气接壤，设置为自由边界；实际情况中坝底与无限大的坝基接壤，因此，设置为无反射边界，以真实地反映爆炸波在边界上的透射情况。为了便于不同计算工况的比较，选取了坝体中的参考单元H1830与H1569进行跟踪，记录该单元的孔压时程曲线和土体密度变化时程曲线等。

图2 有限元模型Fig.2 Finite element model

3.1 考虑炸药埋深变化的弹坑数值模拟

3.1.1 工况介绍

将药包中心距离坝顶的距离定义为炸药埋置深度，在坝顶以上为负，坝顶以下为正。工况分别设置为：①埋深为-0.2 m，即药包位于坝顶表面，为触地爆炸工况；②埋深为0.3 m；③埋深为0.6 m；④埋深为0.8 m；⑤埋深为1.3 m。炸药材料参数如前所述。土体的力学行为通过前述的MAT_FHWA_SOIL材料模型来描述，结合相关经验[13－14,18]，采用的主要参数取值分别为：密度为1.8 g/cm3，土粒相对密度为2.70，体积模量为68.6 MPa，剪切模量为19.4 MPa，内摩擦角为30°，黏聚力为27.3 kPa，屈服面修正系数为34.5 kPa，含水率为0，孔压系数Ksk为0，即在炸药埋深变化的研究中采用干土，不考虑孔隙水压力的发展，进行总应力分析。

3.1.2 结果分析

以坝中轴线为对称轴，将原来的半模型映射为完整模型，不同炸药埋深情况下的堤坝爆炸弹坑形态如图3所示。图中深色部分为坝体形态，中轴线上的小方形表示炸药埋设位置。对模拟结果具体分析如下：

（1）当药包位于坝顶表面时（图 3(a)，触地爆炸），大部分爆炸能量都逸散到空气中去，作用于土体的能量相对较少，因此，只是由于药包下方的土体介质被压实而在坝顶形成一个浅坑。

（2）当药包完全埋置时，由药包下方土体的压实和药包上覆土层的抛掷共同作用产生弹坑。爆炸发生后，压力波从药包中心迅速传出，爆心周围出现一个近球形的空腔向外扩张，药包下方土体形成密实的压缩带。压力波到达坝顶自由面后以稀疏波的形式反射回来，其拉伸作用使得药包上覆土层变得疏松；稀疏波的拉伸以及压力波和爆生气体的推动，导致坝顶不断向上隆起出现鼓包；随着爆炸产物的进一步发展，爆腔两侧的土体开裂扩展到坝顶，与地表连通，爆炸产物便沿着这些裂隙喷发，并携带破碎土块一齐向上方冲出，于是出现了弹坑的雏形（如图 3(b)、3(c)、3(d)所示）；此时，由于炸药埋深的增加，释放到空气中的爆炸能量逐渐减少，而直接作用于土体的能量则逐渐增大，因此，弹坑尺寸随着药包埋置深度的加深而增大。

（3）当药包埋深超过某一深度时，释放到空气中的能量已经忽略不计，埋深的继续增加对其不再产生影响；而此时，由于上覆土体压力较大，主要是由于爆炸的“扩腔”效应使得土体向外扩张，即在爆炸冲击波作用下，土体质点获得速度，沿径向产生位移，炮孔空腔持续扩大；随着冲击波的向外传播，能量迅速衰减，当冲击波到达压缩区界面时，对土体的冲击压缩过程结束，爆腔扩展至极限值，最终在坝顶产生鼓包，没有出现明显的弹坑，而在坝体中则形成一个封闭的地下爆腔（见图3(e)）。

在实际爆炸弹坑试验[19]中，部分土体还会回落，形成一个“乱石井”。本文尚未模拟出土块的回落，但通过多物质ALE法，已经较好地还原了不同埋深条件下爆炸成坑的宏观现象，反映了数值手段的正确性。因此，在进行堤坝抗爆安全设计或爆破拆除堰塞坝时，能够通过数值研究手段预先判断坝体被爆后的形态，为安全设计方案或爆破方案提供一定的指导措施。

在计算中选取了单元H1569（压密区）与单元H1830（空腔区）进行跟踪，记录该单元处的土体密度时程曲线，参考单元的所在位置如图3(d)所示。图4为炸药埋深为0.8 m的工况下，这两个单元中的土体密度变化的对比，由图可以看出，药包下部的单元由于被压实，土体密度变大，而爆腔区中的单元受到爆炸剧烈的冲击压缩作用后，质点向外发生强烈的位移，土体物质被爆炸产物带走，形成空腔，因此，此单元的土体密度在短暂上升之后迅速降为0。

图5为不同炸药埋深工况下堤坝弹坑尺寸的变化关系。随着埋深的增加，弹坑的直径和深度都在加大（见图3(a)～(d)）；但当埋深增加到1.3 m时，无法产生弹坑，而是形成一个封闭的爆腔（见图3(e)）。由于本文采用的是准二维模型，对装药量的影响尚无法进行完全符合实际的定量分析，所以只能够定性地说明不同炸药埋深情况下堤坝爆炸弹坑的变化趋势。

图4 爆腔与压密区土体密度变化(炸药埋深0.8 m)Fig.4 Variations of soil density at explosion cavity and compacting region (burial depth=0.8 m)

图5 炸药埋深与弹坑尺寸关系Fig.5 Relationships of burial depth of explosive and crater size

3.2 考虑土体孔压效应的弹坑数值模拟

3.2.1 工况介绍

几何模型皆为坝顶触地爆炸工况中的模型，即药包中心距离坝顶以上 0.2 m。为了研究爆炸荷载作用下含水土体孔压的上升、有效应力的下降对堤坝弹坑形态的影响，将 MAT_FHWA_SOIL模型中的孔压系数分别设置为 Ksk=0（干土）、1、2、5 MPa，其余参数值不变。如前所述，Ksk越大，材料模型中计算的孔隙水压力就上升得越快，以此来研究土体孔压发展对爆炸成坑效应的影响。

3.2.2 结果分析

图6为不同孔压上升程度下坝顶触地爆炸形成的弹坑形态。由图可以看出，随着孔压系数Ksk的增大，爆炸弹坑在尺寸上也呈加大的趋势，如图 7所示，当孔压系数Ksk=0、1、2、5 MPa时，堤坝产生的弹坑直径分别为1.50、1.60、1.68、3.00 m，弹坑深度分别为0.60、0.70、0.75、0.92 m，说明在爆炸的强烈动荷载作用下，土体的孔压上升导致了土体抗剪强度的下降，加剧了爆炸过程中的土体变形，因此，在爆炸加载下，含水土体较干土会产生更大的弹坑；孔压上升得越剧烈，这种增大就越强烈。

为了直观地显示孔压系数Ksk对土体孔压增长的影响，图 8给出了不同孔压系数Ksk下坝体中参考单元 H1569（见图 3）的超孔压时程曲线，可以看出，当Ksk=0时，超孔压为0，而Ksk越大，单元的超孔隙水压力幅值也越大。

因此，在进行爆破拆除堰塞坝等工程应用时，还必须考虑含水土体的孔压增长对爆破效果产生的影响，可以通过数值手段预先对爆破拆除效果进行分析，合理地控制爆破参数，以避免堰塞坝的坝基或坝体液化而导致的突然溃坝和引发次生灾害[2]。同时可以推断，在爆炸荷载作用下如果堤坝的开裂深度低于水位线，不仅可能导致堤坝失去使用功能，而且将可能造成洪水倾泻乃至溃坝的严重后果，尤其是对于筑坝材料处于饱和状态的土体，在爆炸强烈动荷载作用下产生的超孔隙水压力将进一步加剧爆炸成坑的过程。因此，在遭遇战争、爆炸袭击之前，有必要迅速降低库水水位，以防库水倒灌，引起重大洪灾；同时还可以根据数值手段大致地判断弹坑可能的形状和尺寸，对降水措施提供一定的指导。

图6 不同孔压上升程度下的爆炸弹坑形态Fig.6 Blasting craters under different values of Ksk

图7 不同孔压系数Ksk下的爆炸弹坑尺寸比较Fig.7 Comparison of blasting crater size under different values of Ksk

图8 不同孔压系数Ksk下的土体超孔隙水压力时程比较Fig.8 Comparison of excess pore water pressure time histories under different values of Ksk

4 结 论

（1）利用多物质ALE法的物质流动变形，能够较好地模拟出堤坝在爆炸荷载作用下的弹坑形态。

（2）炸药埋深对堤坝爆炸弹坑形态有重要影响，在埋深较浅时，弹坑的尺寸随着炸药埋深的增加而增大；但当埋深超过某一深度，则无法形成弹坑，而是由于爆炸“扩腔”效应产生一个封闭的爆腔。

（3）爆炸荷载作用下含水土体的孔隙水压力增长对弹坑的形成具有重要影响，孔压上升程度越强，产生的弹坑尺寸越大。

由于本文采用的是准二维模型，计算结果与实际工程的三维工况存在差别，并且尚未考虑抛掷土体回落并覆盖至实际弹坑而导致的可见弹坑与实际弹坑的区别，因此，只能定性地说明爆炸成坑的宏观现象；此外，对堤坝爆炸成坑后的坝体稳定性还未做深入分析。这些都需要进一步开展工作。

[1]梁向前,崔亦昊,魏迎奇. 工程爆破在堰塞湖处理中的应用及实例[J]. 水力发电,2009,35(10): 88－90.LIANG Xiang-qian,CUI Yi-hao,WEI Ying-qi. Treat-ments of landslide dam and its cases by blasting[J].Water Power,2009,35(10): 88－90.

[2]张文煊,吴新霞. 堰塞湖抢险爆破对坝基及保留堰体的液化影响分析[J]. 工程爆破,2008,14(3): 21－23.ZHANG Wen-xuan,WU Xin-xia. Analysis of impact on liquefying of dam foundation and retained weir body in emergency blasting of dammed lake[J]. Engineering Blasting,2008,14(3): 21－23.

[3]SCHMIDT R M,HOLSAPPLE K A. Theory and experiments on centrifuge cratering[J]. Journal of Geophysical Research,1980,85(B1): 235－252.

[4]穆朝民,任辉启,辛凯,等. 变埋深条件下土中爆炸成坑效应[J]. 解放军理工大学学报(自然科学版),2010,11(2): 112－116.MU Chao-min,REN Hui-qi,XIN Kai,et al. Effects of crater formed by explosion in soils[J]. Journal of PLA University of Science and Technology (Natural Science Edition ),2010,11(2): 112－116.

[5]穆朝民,任辉启,李永池,等. 变埋深条件下饱和土爆炸能量耦合系数的试验研究[J]. 岩土力学,2010,31(5):1574－1578.MU Chao-min,REN Hui-qi,LI Yong-chi,et al.Experiment study of explosion energy coupling coefficient with different burial depths in saturated soils[J].Rock and Soil Mechanics,2010,31(5): 1574－1578.

[6]SAUSEVILLE M J. Explosive cratering using a geotechnical centrifuge on dry,saturated,and partially saturated earth embankments and dams[D]. New York:Rensselaer Polytechnic Institute,2005.

[7]ZIMMIE T F,DE A. Centrifuge modeling to study stability of dams[C]//Proceedings of the Annual Conference of the Association of State Dam Safety Officials. Seattle: Association of State Dam Safety Officials,1996: 630－639.

[8]ZIMMIE T F,SAUSVILLE M J,SIMPSON P T,et al.Blasting studies of dams using the geotechnical centrifuge[C]//Proceedings of the ASDSO Annual Conference on Dam Safety. New Orleans: Association of State Dam Safety Officials,2005: 438－445.

[9]刘军,刘汉龙,张正珺. 爆炸荷载下土石坝动力响应特征的数值模拟[J]. 防灾减灾工程学报,2010,30(1): 10－16.LIU Jun,LIU Han-long,ZHANG Zheng-jun. Numerical simulation of dynamic response of an earth and rock-filldam to a blast loading[J]. Journal of Disaster Prevention and Mitigation Engineering,2010,30(1): 10－16.

[10]柳锦春,方秦,还毅,等. 炸药地面接触爆炸下土中感生地冲击的实用计算方法[J]. 解放军理工大学学报(自然科学版),2010,11(2): 121－124.LIU Jin-chun,FANG Qin,HUAN Yi,et al. Practicable calculating method of indirect ground shock in soil at surface contact explosion[J]. Journal of PLA University of Science and Technology (Natural Science Edition ),2010,11(2): 121－124.

[11]FERRERO V H. Further checks for scaling effects on explosion induced craters[D]. Baltimore: University of Maryland,College Park,MD. 1988.

[12]贾祖朋. 基于MOF界面重构的多物质ALE方法[J]. 计算物理,2010,27(3): 353－360.JIA Zu-peng. A multi-material Arbitrary Lagrangian-Eulerian method based on MOF interface reconstruction[J]. Chinese Journal of Computational Physics,2010,27(3): 353－360.

[13]HALLQUIST J Q. LS-DYNA Keyword User’s Manual(971)[M]. California: Livermore Software Technology Corporation,2007.

[14]WAYNE Y Lee. Numerical modeling of blast-induced liquefaction[D]. Provo: Brigham Young University,2006.

[15]KLISINSKI,M. Degradation and Plastic Deformation of Concrete[D]. Warsaw: Polish Academy of Sciences,1985.

[16]KLISINSKI M,MROZ Z. Description of inelastic deformation and degradation of concrete[J]. International Journal of Solids and Structures,1988,24(4):391－416.

[17]LEWIS B A. Manual for LS-DYNA soil material model 147[R]. Washington,D C: Federal Highway Administration,2004.

[18]刘晶波,杜义欣,闫秋实. 地下箱形结构在爆炸冲击荷载作用下的动力反应分析[J]. 解放军理工大学学报(自然科学版),2007,8(5): 520－524.LIU Jing-bo,DU Yi-xin,YAN Qiu-shi. Dynamic response of underground box structures subjected to blast load[J].Journal of PLA University of Science and Technology(Natural Science Edition ),2007,8(5): 520－524.

[19]亨利奇. 爆炸动力学及其应用[M]. 熊建国译. 北京:科学出版社,1987.