远场水下爆炸冲击作用下重力坝动力稳定性的简化时程分析方法

黄宇峰，卢文波，王高辉，陈 明，严 鹏

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

高坝大库由于其显著的政治经济效益，极易成为局部战争[1-2]以及恐怖袭击[3-4]的首要目标，高坝大库一旦失事，上游库区大量洪水下泄，将会给下游的人民带来不可估量的灾难。重力坝作为一种被广泛采用的坝型，其抗爆性能以及安全评价越来越受到重视。

数值模拟和模型试验方法在结构抗爆的研究中得到了广泛的应用[5-10]，国内外学者采用这两种方法，对近场和接触水下爆炸作用下重力坝的动态响应[11-13]、毁伤特性[14-15]以及破坏模式[16-19]进行了大量的研究。胡晶等[12]进行了挡水结构水下爆炸作用的离心模型试验，分别分析了水下爆炸冲击波及气泡脉动对结构的作用机理；王高辉等[15]对水下爆炸冲击作用下重力坝的毁伤机理进行了研究，并据此得到了混凝土重力坝的毁伤模式；LI 等[19]研究了水下接触爆炸下重力坝的破坏模式，探讨了静水压力以及起爆位置对重力坝的损伤特性的影响。水下近场以及接触爆炸工况下冲击波峰值压力大，但衰减得很快，因此重力坝受到的破坏主要为局部毁伤破坏，主要表现为迎爆面及附近的爆炸成坑破坏及冲切破坏、下游坝面的拉伸破坏等，不会对重力坝的整体稳定性造成较大的威胁；而远场起爆工况下，虽然坝面受到的峰值压力较小，但由于其衰减较慢，作用时间更长，影响范围更广，对重力坝坝的整体稳定性造成的威胁更大。因此，需重点关注重力坝在远场水下爆炸冲击下的稳定安全评价问题。

近年来有关爆炸冲击等荷载作用下的结构动力稳定性能的研究取得了较大的进展[20-22]，在重力坝的稳定计算方面，目前的研究主要集中于复杂地质条件下的深层抗滑稳定计算以及地震作用下重力坝的动力稳定性计算，其主要计算方法有刚体极限平衡法以及有限元法。刚体极限平衡法[23]将大坝视为刚体，以滑移体上的抗滑力与滑动力之间的关系来评价重力坝的稳定性，其力学概念明确，计算方法简便，是重力坝稳定分析的主要方法。近年来有学者提出了一种动力有限元时程法，以有限元方法计算区域内单元的应力，积分得到滑移面上的滑动力与抗滑力，再结合极限平衡理论，计算得到稳定安全系数，该方法在重力坝抗震动力稳定分析中应用比较广泛[24-25]。

综上，目前尚缺乏对远场水下爆炸冲击作用下重力坝稳定性的研究。本文正是聚焦于此问题，综合考虑了远场水下爆炸冲击波的传播及其与重力坝的相互作用过程，基于波动理论以及平面波假设，建立了远场水下爆炸冲击荷载作用下重力坝动应力场的简化理论计算模型，并结合时程分析方法，建议了一种远场水下爆炸冲击作用下重力坝沿建基面的稳定安全系数的计算方法，为重力坝的抗爆安全评价提供理论参考。

1 远场水下爆炸冲击波的传播及其与重力坝的相互作用过程

本文研究中的远场起爆工况为起爆距离大于50 倍的药包半径的水下起爆工况。炸药在上游库区起爆后，会以球面波的形式向四周传播。当冲击波传播至上游坝面后，会与坝面相互作用，一部分冲击波会透射进入坝体内部，以应力波的形式继续向下游坝面传播； 当应力波传播至下游坝面时，会在下游坝面发生反射，形成反射纵波与反射横波；反射波继续向坝上游侧传播，当其传播至上游坝面或是建基面等结构面时，又会发生透反射，如此循环往复，直至波的能量完全衰减。远场水下爆炸的冲击波的传播及其与坝体的相互作用过程十分复杂，为了便于研究，本文对其进行了简化：假设冲击波传播至上游坝面时，以平面波的形式作用于坝体[26]；由于波在传播过程中不断的衰减，特别是在坝体混凝土与其他介质分界面处发生透反射后，波的能量将大幅削弱，因此在研究时只考虑冲击波在坝面的透反射以及随后应力波在下游坝面的反射对坝体稳定性的影响，忽略反射波在坝基截面等其他分界面的透反射以及之后坝内波的传播等因素。

简化之后的远场水下爆炸冲击波的传播及其与重力坝的相互作用过程如图1 所示。

图1 远场水下爆炸冲击波的传播及其与重力坝的相互作用过程示意图Fig.1 Schematic diagram of far-field underwater explosion shock wave propagation and its interaction with gravity dam

1.1 水下爆炸冲击波与坝面的相互作用

自由场水下爆炸冲击波峰值压力与时程曲线一般采用Cole 经验公式计算[27]：

传播至上游坝面的冲击波会在坝面发生透反射，反射压缩波与入射压缩波叠加，导致坝面压力陡增，如图2 所示。根据平面弹性波在不同介质间的透反射公式[28]，大坝表面实际受到的压力为：

图2 爆炸冲击波与坝体的相互作用Fig.2 Interaction between explosion shock wave and gravity dam

式中：pm、分别为大坝表面实际峰值压力以及入射峰值压力；p(t)为大坝表面实际压力时程曲线；ρ1、ρ2分别为水和坝体混凝土的密度；CP1、CP2分别为波在水中和混凝土中的波速。波在水中的波速可取1500 m/s，在混凝土中的波速CP根据式(6)计算得到，E、ρ、µ分别为坝体混凝土的弹性模量、密度和泊松比。

1.2 爆炸应力波在下游坝面的反射

透射进入坝体内部的应力波为纵波，当其传播至下游坝面处时，透射波作为新的入射纵波IP，在下游坝面处发生反射，形成反射纵波RP和反射横波RS，如图2 所示。

由于波传播过程中能量的耗散，传播至下游坝面的应力波可视为弹性平面波。为了得到应力波在下游坝面处的反射系数，基于弹性平面波的假设，本文对入射纵波IP、反射纵波RP和反射横波RS在下游坝面处引起的单元体应力状态进行了分析，如图3 所示。图中σx1、σx2、σx3分别为入射纵波、反射纵波以及反射横波引起的正应力，τxy1、τxy2、τxy3分别为入射纵波、反射纵波以及反射横波引起的剪应力。

图3 不同应力波作用下交界面处单元体应力状态Fig.3 Stress state of unit body at the interface under the action of different stress waves

远场起爆条件下爆炸荷载周期较长，一般可达十几甚至几十毫秒，其波长与重力坝的截面宽度较为接近，其引起的结构微元的惯性效应较小，为了简化计算，本文在分析不同应力波作用下交界面处单元体的应力状态时，忽略了惯性效应的影响。

入射纵波所引起的三角形单元应力状态如图3(a)所示，根据静力平衡条件[29]可得：

式中：α 为入射角；AC=ABcosα，BC=ABsinα。

将AC、BC代入式(7)，解得：

对单元体ADE，根据平衡方程可得反射纵波引起的正应力与剪应力分别为：

对单元体AFG，根据平衡方程可得反射横波引起的正应力与剪应力分别为：

式中：α′为纵波反射角；β 为横波反射角，根据波传播的Snell 定律，入射角与反射角之间满足以下关系：

式中，CS为反射横波波速，其计算式为：

不考虑下游库水，将重力坝下游坝面视为自由面，在自由面上有 σx=0、τxy=0，即：

将式(8)～式(10)代入式(13)，即可得到入射纵波在下游坝面的纵波反射系数k1以及横波的反射系数k2：

式中，A1、B1、C1、D1各系数的表达式为：

2 坝内动应力计算

应力波在坝体内传播的过程中会发生能量的衰减，其应力峰值会随着传播距离的增大不断减小，由于岩石和混凝土的物理性质较为接近，可参考爆炸应力波在岩石介质中峰值压力的衰减公式：

式中：R为起爆距离，即爆源到重力坝上游坝面的距离；pm为坝面压力峰值；x为应力波在坝内传播的距离；γ 为坝体混凝土中峰值压力的衰减系数，在弹性区可取1[30]。

根据前文的假设可知，对于坝体内任意一点M，该处的动应力应为从上游坝面透射进入坝体内部的P 波以及下游坝面反射形成的P 波和SV 波三者的叠加，如图4 所示。

图4 坝体动应力计算示意图Fig.4 Schematic diagram of dam dynamic stress calculation

入射纵波IP在M点引起的动应力为：

式中：LOM为入射纵波传播路径OM的长度；t1为入射纵波从坝面传播至M点所需要的时间，计算式为：

从O1 处入射的纵波，在下游坝面J 处发生反射后，形成的反射纵波RP，经JM 传播至M 处，在反射纵波的作用下，M 点的动应力为：

从O2处入射的纵波，在下游坝面K处发生反射后，形成的反射横波RS，经KM传播至M处，在反射横波的作用下，M点的动应力为：

在M点的位置给定以后，根据几何关系，便可确定入射纵波、反射纵波以及反射横波的传播路径，其表达式为：

式中：θ 为混凝土重力坝下游坝坡坡角；ε 为反射纵波与水平面之间的夹角，；φ为反射横波与水平面之间的夹角，

在分别分析得到入射P 波、反射P 波以及反射SV 波在M点引起的动应力后，将其分别向x方向以及y方向分解，得到爆炸应力波作用下坝体质点的正应力σy与剪应力τyx。

对于入射纵波IP，其表达式为：

对于反射纵波RP，其表达式为：

对于反射横波RS，其表达式为：

将三组波所引起的应力分别叠加，即可得到应力波作用下M点的正应力与剪应力的表达式：

3 重力坝远场爆炸动力稳定性简化时程分析方法

根据第1 节和第2 节推导得出的远场水下爆炸冲击作用下重力坝动应力场计算公式，即可求解重力坝建基面上的动应力场，结合时程法，就能得到重力坝在任意时刻沿建基面的抗滑稳定安全系数。

混凝土重力坝远场爆炸动力稳定性简化时程分析方法的具体步骤如图5(a)～图5(d)所示。

图5 重力坝远场爆炸动力稳定性简化时程法计算流程图Fig.5 Calculation flow chart of simplified time history method for far-field explosion dynamic stability of gravity dam

1) 首先根据给定的起爆药量、起爆距离等参数，计算得到水下爆炸冲击荷载，考虑波在坝面的透反射作用，计算得到坝面冲击荷载时程曲线。

2) 计算应力波在下游坝面处的反射系数，得到透射波σIP、反射纵波σRP以及反射横波τRS的应力时程表达式。

3) 分别计算入射纵波、反射纵波和反射横波在建基面上引起的动应力，然后将三者引起的动应力叠加，得到建基面上任意质点的正应力σy与剪应力τyx。要注意的是，上述的正应力与剪应力是由水下爆炸冲击荷载引起的，除此之外，建基面上的总应力还应该考虑静力荷载作用下(重力G、静水压力荷载pw以及扬压力荷载pu)产生的静应力和。建基面上的静应力可根据材料力学法确定。这样，就能得到建基面上总的应力分布。

4) 将建基面分为n个条块，每个条块的长度记为li，取每个条块中点处的应力代表整个条块的应力，于是对于每个条块，其滑动力为抗滑力为，再将每一份块体的抗滑力与滑动力累加，即可得到沿建基面整体的滑动力与抗滑力，两者的比值即为重力坝沿坝基面的抗滑稳定安全系数：

式中，f′、c′分别为坝体混凝土与坝基接触面的抗剪断摩擦系数及粘聚力。

4 混凝土重力坝动力稳定性分析

4.1 计算模型与参数

本文选取某混凝土重力坝典型坝段作为动力稳定性分析对象，该重力坝高度为80 m，坝基底部宽度为66 m，下游坝面折坡比为1∶0.8，上游正常蓄水位深度为75 m，坝段宽度为10 m。为验证本文所建议方法的可靠性，建立了混凝土重力坝数值计算模型，如图6 所示。本文采用声固耦合方法，声单元模拟水下爆炸冲击波在库水中的传播，实体单元模拟大坝、地基。坝体混凝土、基岩均假定为线弹性材料，弹性模量为30 GPa，泊松比为0.167，密度为2400 kg/m3，单位类型为C3D8R；库水选用声学单元AC3D8R，其体积模量为2.14 GPa，密度为1000 kg/m3。无反射边界条件是通过acoustic impedance(声阻抗)的方式施加于水网格的边界上，因为水体边界是平面的，所以施加平面无反射边界。本文荷载施加采用输入压力时程曲线的方式，入射波荷载以平面波的形式施加于库坝耦合面，入射波压力时程曲线根据式(1)～式(3)计算所得，如图7 所示。

图6 混凝土重力坝数值计算模型Fig.6 Numerical calculation model of concrete gravity dam

图7 上游坝面入射压力时程曲线Fig.7 Time history curve of incident wave pressure at the upstream dam surface

本文在进行稳定计算时，考虑的荷载组合为：上游坝面静水压力、扬压力、重力坝自重以及水下爆炸冲击荷载；坝体混凝土与坝基接触面的抗剪断摩擦系数f′取0.9、粘聚力c′取0.8 MPa。本文选取的计算工况为起爆药量W为5000 kgTNT，起爆距离R为100 m。

4.2 方法可靠性验证

为了验证本文所建议的简化时程分析方法的可靠性，基于ABAQUS 数值模拟软件对远场水下爆炸作用下重力坝内部的应力波传播过程进行了模拟，得到了重力坝建基面质点动应力以及沿建基面的稳定安全系数的数值解，并与弹性平面波假设下的解析解进行了对比。由于数值计算中坝体以及基岩材料选取的是线弹性模型，在数值计算中无法模拟出应力波在传播过程的能量损耗带来的应力波峰值的衰减，因此选取衰减系数γ 为0 时的解析解与数值解进行对比，如图8 所示。

图8 解析解与数值解的对比Fig.8 Comparison between analytical solution and numerical solution

根据图8(a)、图8(b)可以发现，解析解与数值解在前期偏差较小，在后期数值解的波形较乱，与解析解之间的偏差较大。根据图8(c)可以发现，稳定安全系数的数值解与解析解在30 ms之前吻合得较好，两种方法计算得到的稳定安全系数的变化趋势一致，在数值上的偏差也较小，前20 ms 两者的偏差在20%以内，20 ms～30 ms 的偏差在30%以内。在30 ms 后，数值解出现了与解析解不一样的变化规律，这是因为此时上游坝面反射的应力波传播至坝基截面，在数值计算中，由于没有考虑波在传播过程中的衰减，坝基截面上产生了较大的剪应力，所以导致稳定安全系数显著下降。但是在实际情况中，波在传播过程中存在着能量的耗散，上游坝面反射形成的应力波在传播至坝基截面后压力很小，对重力坝动力稳定性的影响较小，因此在解析解中忽略了上游坝面反射波的影响。这就导致了在后期数值解与解析解之间的差距较大。总体看来，解析解与数值解在变化趋势上具有良好的一致性，在数值上大部分时间两者的偏差较小，从侧面说明了该解析方法具有较好的可靠性。

4.3 大坝动力稳定性分析

图9 给出了衰减系数γ 取1 时重力坝建基面形心处的动应力时程曲线。从图9 可以看出，坝基截面形心处质点的动应力时程曲线可以明显分为以下3 个阶段：Ⅰ：透射纵波作用下的受压阶段；Ⅱ：反射纵波作用下的拉剪阶段；Ⅲ：反射横波作用下的压剪阶段。透射进入坝体内部的纵波为压缩波，在其作用下，质点将处于受压阶段；当应力波传播至大坝下游面后，会发生反射，形成一组纵波和一组横波，纵波为拉伸波、横波为剪切波。反射纵波传播速度快，因此会先于反射横波传播至建基面形心处，在其作用下，该处质点会产生拉应力以及剪应力，该阶段质点处于拉剪阶段；之后，反射横波也传播至该处，导致该处质点产生较大的剪应力与压应力，该阶段质点处于压剪阶段。

图9 建基面形心处质点动应力时程曲线 (γ=1)Fig.9 Time history curve of particle dynamic stress at the centroid of foundation surface (γ=1)

图10 给出了衰减系数取1 时建基面上抗滑力、滑动力以及稳定安全系数随时间变化的规律。根据图10 可知，重力坝的稳定安全系数变化呈现出先增大后减小再增大并最终保持恒定的变化趋势。初始时刻重力坝只受到透射压缩波的作用，透射压缩波在坝基截面引起压应力，使得建基面上的抗滑力增大，从而导致初始时刻重力坝的稳定安全系数会略大于静力稳定安全系数。但反射形成的拉伸波传播至建基面后，会导致建基面产生拉应力与剪应力，而当剪切波传至建基面后，又将会引起较大的剪应力以及压应力，在两组反射波的共同作用下，建基面上的滑动力显著增加，同时抗滑力也显著下降，重力坝的稳定安全系数不断减小，在22 ms 时达到最小值。之后，随着应力波的衰减，大坝的稳定安全系数不断增大，逐渐趋于静力稳定安全系数。

图10 稳定安全系数时程曲线图(γ=1)Fig.10 Time history curve of stability safety factor (γ=1)

5 结论与展望

本文基于平面波假设和波动理论，对远场水下爆炸冲击作用下重力坝的稳定安全问题进行了研究，得到以下结论：

(1) 基于远场水下爆炸冲击波的传播及其与坝体的相互作用过程的分析，建立了重力坝爆炸动应力场的简化理论计算模型，最终建议了一种远场水下爆炸冲击作用下重力坝动力稳定性的简易时程分析方法，并得到了动力有限元数值方法的验证。该方法计算简便，能较为快速地计算得到重力坝沿建基面的动力稳定安全系数时程曲线，具有一定的可靠性。

(2) 在远场水下爆炸冲击作用下，重力坝沿建基面的抗滑稳定安全系数呈现出先增大后减小再增大并最终保持恒定的变化趋势。反射横波与反射纵波的共同作用会导致建基面的滑动力增大，抗滑力下降，使得重力坝的稳定安全系数显著减小。

本文仅仅只是对水下爆炸打击下重力坝动力稳定性问题进行了一个初步的研究，在实际的远场水下爆炸打击下，重力坝坝基会在冲击荷载的作用下开裂，基岩的力学性质也会在应力波扰动下发生变化，如何在稳定分析中将这些因素加以考虑，仍有待进一步的研究。