水下爆炸荷载作用下圆柱结构反射压力解析计算方法研究

闫秋实，张志杰，王丕光，宝 鑫，李述涛

(1. 北京工业大学城市与工程安全减灾教育部重点实验室，北京 100124；2. 清华大学土木工程系，北京 100084)

随着交通需求的迅速提高，渡海桥梁、高桩码头等近海结构大量兴建。然而现今恐怖事件频发，近海结构因其建造复杂，易炸难修，常是恐怖袭击或军事打击的重点目标，其存在遭受水下爆炸破坏的可能性，因此，开展水下爆炸荷载的研究对于结构进行抗爆防护、毁伤预测和控制具有重要意义[1]。

水下爆炸的系统研究最早可追溯Cole[2]在1948年出版的《Underwater Explosion》专著。在书中Cole[2]基于大量水下爆炸试验，最早系统总结了水下爆炸研究成果，对水下爆炸基本规律，现象和试验研究方法做了详细的叙述，并基于基尔克乌特-别泽理论提出了自由场的半经验半理论压力计算公式。Bobrovskii等[3]在Cole[2]的研究基础上进行补充，重点分析了水面、水底和结构物对冲击波产生的反射、绕射和空化效应，并对Cole[2]压力计算式进行完善。目前水下爆炸载荷的研究方法主要可分为试验方法、解析方法[4]和数值方法，数值方法又可分为双渐进近似法(DAA)[5]、有限元法(FEM)、有限差分法(FDM)[6−7]、光滑粒子法(SPH)[8]和间断伽辽金法(DG)[9]等。基于这些研究方法，国内外研究人员得到水下爆炸载荷的传播变化规律。Geers和Hunter[10]基于DAA法提出了一个考虑气泡震荡上浮的压力数值模型，该模型被广泛应用于数值计算。伍俊等[11]和明付仁[12]通过数值模拟对Cole[2]提出的峰值压力预测计算式进行了修正。

水下环境复杂，冲击波在传播过程中的规律变化受多种因素影响。王高辉等[13]和明付仁等[14]分析了水下爆炸冲击波在近边界面上的传播特性，研究表明：当冲击波在传播过程中与结构发生接触时，由于水和结构的波阻抗不同，结构表面会反射产生稀疏波，入射波与反射波在结构表面叠加促使压力出现骤然上升。Zhuang等[15−16]进行了一系列水下爆炸单桩试验，依据试验结果得到了冲击波压力沿桩身的分布规律，并提出便于工程应用的压力工程算法。Wang等[17]基于水体为可压缩、无旋和无粘流体的假设，采用时域人工边界方法提出了远场水下爆炸冲击波下圆柱表面反射压力的子结构分析模型。Iakovlev[18]基于拉普拉斯变化和分离变量法推导了二维壳体受冲击波作用的散射和绕射波场。目前，对水下爆炸荷载作用下结构反射压力的研究已取得一定成果，但多数采用数值或试验方法，这两种方法尽管计算结果的精度高，但费时低效不利于快速应用计算。鲜有见到采用解析方法处理作用于近海结构中常用的圆柱结构上的反射压力。尽管解析方法仅适用于圆柱等规则结构，但其计算效率高，在一定程度上有助于快速、大量计算结构受冲击波作用下受损和破坏情况，这对其进行有效抗爆防护，快速修复具有一定指导意义。此外，水体压缩性也是影响水下爆炸荷载作用下圆柱反射压力的重要因素，目前研究还缺乏对其相关进行定量分析。

本文拟基于绕射波浪理论，在忽略水体压缩性的假设条件下，根据水体控制方程和边界条件，在柱坐标系下通过分离变量推导水下爆炸冲击波作用下圆柱表面分布反射压力的时域解，并根据Wang等[17]的时域人工边界方法，对比考虑水体压缩性和忽略水体压缩性的总动水力峰值差异，分析水体压缩性对反射压力的影响。

1 水下爆炸荷载作用下圆柱结构的反射压力解析解

当炸药在水中发生爆炸时，产生的爆炸冲击波会以球面波的形式在水中传播并与水中结构发生相互作用。水中结构会对冲击波传播路径产生影响，出现冲击波的反射和绕射现象，而水中结构也会受水下爆炸冲击波的作用发生振动扰动周围水域产生辐射波，如图1所示。因此，水下爆炸冲击波对圆柱结构作用可考虑为入射波、反射波和结构受荷振动产生的辐射波的叠加作用。本文假设齐水面圆柱为刚性，首先基于水下爆炸冲击波自由场压力经验公式，为圆柱结构施加入射压力，接着基于绕射波浪理论推导柱坐标系下反射压力解析解。

图1 冲击波与圆柱结构产生相互作用Fig. 1 The interaction between shock wave and circular cylinder

1.1 水下爆炸荷载作用下冲击波入射压力

水下爆炸荷载作用下齐水面等截面圆柱与水相互作用分析模型如图2所示，图2中：a为圆柱半径；h为水深(柱高)。柱坐标系和直角坐标系的z轴沿柱中心向上，坐标系原点位于柱底部中心处。水的密度和波速分别取为1000 kg/m3和1438 m/s。假定柱底地基为刚性，不考虑冲击波受地基和水面反射的影响。爆距与炸药半径的比值为比例爆距Dr[19]，当10＜Dr≤25时，水下爆炸属于近场爆炸，为非接触爆炸，主要受冲击波作用；当Dr＞25时，水下爆炸属于远场爆炸。

图2 等截面圆柱与水相互作用示意图Fig. 2 The system of vertical cylinder interaction with fluid

当炸药在水下发生爆炸时，爆炸冲击波以球面波形式向外扩散，Cole[2]将冲击波压力在自由场中的衰减过程表述如下：

式(1)～式(3)中：pi(R,t)为爆炸冲击波压力时程；Pm为爆炸冲击波峰值压力；θ1为时间衰减常数；t0为冲击波到达的时刻；W为炸药药量；K1、a1、K2、a2为对应于TNT炸药的峰值压力和时间衰减常数的经验常数；R为起爆点距某点处的距离，起爆点的坐标可表示为(x0, y0, z0)，则在柱坐标系下R可表示为式(4)的形式。其中，R与t0的表达式如下：

式中，c为水中声速。

1.2 水下爆炸荷载作用下冲击波反射压力

假设水体为不可压缩，则在柱坐标系下不可压缩水体中反射压力ps可通过拉普拉斯方程表示为：

而反射压力在柱坐标系中满足以下边界条件：

1)刚性地面条件：

2)自由表面条件：

3)水体与结构交界面处边界条件：

4)无穷远处辐射边界条件：

通过分离变量法反射压力可表示为：

并将式(11)代入式(6)中则可拆分成如下方程：

式中，r0=rλ。

则式(12)的解可表示为：

式中，Am和Bm表示待定系数。

由式(13)通过边界条件式(7)和式(8)中可有：

式中，λ=(2 j−1)π/2h，且Cj为待定系数。

将无穷远处的辐射边界条件式(10)代入式(14)中则有：

式中，Dm为待定系数。

将上述式(15)～式(17)代入拉普拉斯式(6)中，则反射压力可表示为：

式中，Km(·)为第二类修正的贝塞尔函数。

根据水与结构交界面的边界条件式(9)，式(18)中的Ejm(t)和Fjm(t)均表示待定系数，可表示为：

则总压力p可表示为：

2 方法验证对比及分析

2.1 解析验证

Wang等[20]基于水体控制方程和边界条件，提出了一种应用于计算地震作用下刚性圆柱动水力的高精度高效数值计算模型。本文应用该模型验证提出的反射波解析解，其中验证模型设置为齐水面等截面圆柱，柱高为4 m，半径为0.4 m。炸药药量10 kg，起爆点位于距圆柱中心2.5 m处。

图3表示在2 ms时刻，圆柱中心环向位置上的反射压力对比，由图3可看出解析解与数值解吻合程度较好。

图3 柱中位置反射压力对比Fig. 3 Comparison of scattered wave pressure at column middle position

2.2 对比有限元计算

开展水下爆炸研究常选用有限差分法(FDM)，该方法通常有两种思路。一种是基于声学介质理论，不考虑炸药的爆轰过程，将冲击波直接作用到结构上，通过声固耦合进行动力响应计算；另一种考虑炸药的爆轰过程，建立炸药起爆，可用于模拟冲击波传播和与结构发生相互作用的过程，常用软件有Autodyn、Ls-dyna[21]等。闫秋实等[22]、Yang等[23]和徐强等[24]通过Autodyn分别模拟了水下爆炸荷载作用下混凝土重力坝和钢筋混凝土柱破坏损伤情况，分析了其动力响应和破坏机理。

2.2.1 对比模型设置

本文拟选用Autodyn进行解析计算与有限元计算进行对比。对比模型同样设置为齐水面圆柱，柱高为4 m，直径为0.8 m，网格大小为50 mm，柱底部固定。考虑到模型计算精度要求，水域大小设置为4.5 m×2 m×4 m，并在表面设置流出边界，对流体与结构交界处的流体进行局部加密，网格大小设置为25 mm，对比模型示意如图4所示。炸药选用10 kg的TNT，炸药半径则为0.113 m，起爆点距离结构表面2.5 m，离水面2 m，根据比例爆距计算：2.5/0.113≈22，为非接触爆炸情况，结构主要受冲击波作用。通过Autodyn内置的映射技术，将水下爆炸产生的一维冲击波自由场压力计算结果加载到三维圆柱结构表面，进行压力计算。

图4 有限元对比模型 /mFig. 4 Finite element comparison model

沿柱高每沿0.5 m在迎爆面和背爆面上设置测点，且在柱中位置处沿环向每隔30°设置测点，迎爆面上编号依次为Y1～Y7，背爆面上编号依次为B1～B7，柱中心环向位置编号依次为H1～H10，圆柱测点编号如图5所示。分别比较解析和有限元计算得到的压力时程、冲量时程和峰值压力。

图5 测点布置图Fig. 5 Measuring point layout

2.2.2 参数对比

图6和图7显示了圆柱迎爆面和背爆面上沿柱高度方向分布测点处的峰值压力对比。

图6 迎爆面峰值压力对比Fig. 6 Contrast of peak pressure of cylindrical upset surface

图7 背爆面峰值压力对比Fig. 7 Contrast of peak pressure of cylindrical back explosion

从图6和图7可明显看出，解析和有限元计算得到的反射压力分布规律相同。在迎爆面上，柱的中心位置处(Y4)峰值压力有限元解高于解析解，偏差约为19.61%，且底部处二者误差最大，达到35.6%，造成这些误差是由于解析方法尚未考虑水体压缩性，而冲击波在频谱成分中高频波的占比多，高频波会受水的压缩性影响，而Autodyn计算可考虑水体压缩性并采用显式算法对动力学方程求解，因此，解析解在峰值压力上计算会低于有限元解，并且柱底部存在多次反射波作用，而解析解尚无法考虑刚性地面的反射波多次作用，因此造成底部误差较大；在背爆面上，峰值压力沿柱高度各测点偏差不一，除顶点外平均误差在15.58%左右，柱顶部偏差最大达到38.13%，解析解高估柱顶处的峰值压力。图8表示柱中环向测点的峰值压力对比，从图8中可发现，峰值压力分布规律相同，总体上峰值压力沿环向解析解高于解析解，在Y4等局部测点处解析解在峰值压力上低于有限元解。造成环向上的峰值压力差异是由于本文解析方法基于声学理论，与有限差分算法不同，在环向部分测点上高估峰值压力，从图6和图7中也可看出峰值压力差异分布规律不一。

图8 环向峰值压力对比Fig. 8 Cylindrical hoop peak pressure comparison

图9显示Y4测点处的压力和冲量时程曲线对比图，由图9可知解析方法得到的峰值压力低于有限元计算的峰值压力，误差约为19.61%，对比冲量可发现，解析解相对于有限元解在压力上衰减更快，因此在2.5 ms前冲量达到峰值；而有限元计算的冲量值由于考虑了刚性地面的反射波作用，其值略高于解析值，误差约为6.13%，说明解析解尽管在峰值压力上由于忽略水的压缩性的原因计算结果偏小，但在作用的冲量上与有限元计算相差较小，对结构的作用效果大致相同，因此说明该解析方法可应用于计算，误差在工程允许误差范围内。

图9 Y4测点压力和冲量时程曲线Fig. 9 Y4 measuring point pressure and impulse time history curve

2.3 分布规律分析

在前反射压力时域解的前提下开展了爆炸距离的工况分析，探讨在相同炸药药量情况下，在不同爆距的情况下迎爆面上反射压力的变化规律。

分析模型设置为半径为0.4 m，高为4 m的圆柱，并设置10 kg TNT炸药位于水下2 m起爆，爆距分别设置为1 m、1.5 m、2 m、2.5 m、3 m、3.5 m、4 m和4.5 m，则比例爆距为8.85、13.27、17.7、22.12、26.55、30.97、35.4和39.82。计算结果如下图10和图11所示。

图10显示不同比例爆距下峰值压力沿柱高分布，从图10中可看出比例爆距越小则峰值压力分布越不均匀，柱中位置处峰值压力最大，逐渐向两端减小。随着比例爆距增大，沿柱高分布的峰值压力趋向均匀，当比例爆距大于26.55时，可将冲击波视为平面波，峰值压力沿柱高相同。图11表示反射峰值压力与入射峰值压力比值在柱中和柱底位置处随比例爆距变化情况，从图11中可看出两处典型位置处反射峰值压力与入射峰值压力比值随比例爆距增大而增大。对于柱中位置，当比例爆距大于17.5时，反射峰值压力可视为1.45倍入射峰值压力；对于柱底位置，当比例爆距大于30时，反射峰值压力可视为1.8倍入射峰值压力。

图10 不同比例爆距下峰值压力分布Fig. 10 Peak pressure distribution under different proportions of burst distance

图11 不同比例爆距下典型位置峰值压力比值分布Fig. 11 Distribution of peak pressure ratio under different proportions of burst distance

3 水体压缩性影响讨论

本文中由于暂未考虑水体压缩性的影响，假设水体为不可压缩，而实际情况中水体压缩性会影响波的传播过程。Liaw和Chopra[25]最早研究了水体压缩性对齐水面圆柱地震动水压力的影响，其研究结果表明：对于细长的圆柱可忽略水体压缩性的影响。Tanaka等[26]和Anthony[27]对圆柱动水力研究表明：当考虑水体压缩性时，会存在辐射阻尼效应。而由于辐射阻尼的存在，即使没有结构阻尼，结构达到共振时其动力反应也不会无限增大。因此，需对水下爆炸荷载作用下圆柱反射压力受水体压缩性的影响进行量化分析。

Wang等[17]假设水体为可压缩、无旋和无粘流体，采用时域人工边界方法建立了远场水下爆炸冲击波下圆柱表面反射压力的子结构分析模型。因此，本文基于Wang等[17]的研究方法，计算冲击波作用下的可压缩水体中圆柱表面反射压力，对比分析不同尺寸的圆柱结构宽深比、水深情况下水体压缩性对结构反射压力的影响。

3.1 分析模型设置

现将圆柱宽深比2a/h、水深h和炸药药量W设为变量参数，分别设置圆柱宽深比2a/h为0.1～0.5，水深为5 m、10 m和15 m，以及炸药为10 kg和50 kg，共计60个工况来探究分析，如表1所示。柱高与水面平齐，柱底端固定。TNT炸药起爆点均设置位于柱高1/2，距圆柱表面水平5 m位置处。

表1 工况设置Table 1 Condition setting

将各点处不同时刻的冲击波反射压力通过环向积分和纵向积分，可得到水下爆炸载荷作用下反射波引起的总动水力，通过对总动水力F这一结构全局指标变化来分析判断水体压缩性对于反射压力的影响，积分形式可表示为：

式中：ps为任意时刻圆柱表面某点处反射压力；a为圆柱半径；Ai为高度为第i层处通过环向积分得到某一时刻的反射压力之和；F为某一时刻通过纵向积分得到的水下爆炸载荷作用下反射波引起的总动水力。

3.2 宽深比和水深对反射压力的影响

图12显示了圆柱宽深比和水深对动水力的影响。由图12可看出，圆柱宽深比越大则水体压缩性影响越显著。当圆柱宽深比为0.1较小时，不同水深情况下可压缩水体和不可压缩水体中的动水力时程曲线均较为吻合，动水力的峰值有所差异，说明在圆柱宽深比较小时，动水力受水体压缩性的影响较小；当圆柱宽深比由0.1增大至0.5时可发现，随水深从5 m变化至15 m，动水力峰值差值分别为25%、28.57%和36.36%，说明动水力峰值受水体压缩性影响随水深的增加而增大。当水深越大时，不可压缩水体与可压缩水体的动水力峰值差值越大，动水力时程曲线差异也越大。

图12 宽深比和水深对动水力的影响Fig. 12 The influence of cylinder width to depth ratio and water depth on hydrodynamic pressure

图13为不同工况下考虑水体压缩性和忽略水体压缩性对动水力F峰值影响的比较。图13中纵轴中的F1为不考虑水体压缩性时的动水力峰值，F2为考虑水体压缩性时的动水力峰值。当F1/F2的比值越接近1时，就表明两种水体差异越小，动水力受水体压缩性的影响越小。由图13中可以看出：随水深和宽深比的增大，F1/F2的比值与1差值越大，水体压缩性对反射压力影响较大，说明忽略水体压缩性会使反射压力计算结果偏大。当圆柱的宽深比小于0.2时，在水深为5 m、10 m和15 m情况下，F1/F2的比值均小于1.2，说明对于柔细圆柱结构水体压缩性对结构反射压力的影响较小；当宽深比大于0.2时，在水深为5 m、10 m和15 m情况下，F1/F2的比值基本大于1.2，说明此时水体压缩性对反射压力影响较大，由于水体中辐射阻尼效应，冲击波在传播时耗散了能量。

图13 水体压缩性对反射压力的影响Fig. 13 Water compressibility effects on maximum hydrodynamic force of the elastic cylinders in time domain under different explosive charge

4 结论与展望

本文基于绕射波浪理论，在柱坐标系下通过分离变量法和边界条件根据水体控制方程推导了水下爆炸荷载作用下圆柱反射压力的解析解，并以数值模拟为基准，比较验证了解析解的有效性，比较结果显示解析算法与数值模拟误差在工程允许范围内。其次，探讨了水体压缩性对反射压力的影响，分析表明：当圆柱宽深比小于0.2时，考虑水体压缩性和忽略水体压缩性时对于动水力峰值影响不大，说明对于细长结构的圆柱结构，水体压缩性对于冲击波反射压力影响较小，可按不可压缩水体进行计算。当宽深比大于0.2时，圆柱为矮粗结构，动水力峰值受水体压缩性影响较大。忽略水体压缩性时的动水力峰值显著大于考虑水体压缩性时的动水力峰值，这是由于可压缩水体考虑了水体辐射阻尼效应的影响，冲击波在水中传播发生了能量耗散。因此，当忽略了水体压缩性的时，会造成反射压力计算结果偏大。

本文提出的解析方法在工程误差允许范围内能有效计算圆柱的反射压力，但在峰值压力上受水的压缩性影响还与有限元计算存在差距，此外还尚无法考虑水面稀疏波和刚性地面反射波对结构入射的影响，因此，需进一步优化解析方法，来考虑上述因素的影响达到更精确的计算。