水下爆炸圆柱壳冲击波绕射特性分析

2020-02-08姚熊亮屈子悦张乐王志凯

哈尔滨工程大学学报 2020年1期

姚熊亮，屈子悦，张乐，王志凯

(哈尔滨工程大学船舶工程学院，黑龙江哈尔滨 150001)

在水下爆炸作用下，由于高能炸药的爆炸、冲击波的传播以及复杂的流体-结构相互作用现象的存在，使得结构周围的水动力场变得十分复杂。潜艇作为现代海战的重要组成部分，由于其圆柱形结构的特殊性，在分析水下爆炸冲击波波动效应时，除了要考虑反射、透射冲击波外，在圆柱壳背爆面出现的绕射冲击波也不容忽视。

姚熊亮等[1]使用ABAQUS有限元软件，利用声固耦合法对水下爆炸载荷作用下水面舰艇和潜艇结构的动态响应进行了数值模拟。朱正洋[2]对近地空爆冲击波进行了数值模拟，得出了冲击波在空气介质中的传播过程和作用在结构上的绕射载荷效应。Ofengeim等[3]通过数值方法研究了平面激波在刚性圆柱体上的绕射。年鑫哲等[4]分析了空气冲击波作用于柔性防爆墙的透射和绕射效应，比较了压力波形的变化，得到了墙后压力场的变化规律。邵宗战[5]为研究水下爆炸对鱼雷、水雷等水下目标的破坏机理开展了一次水下爆炸试验，测量了圆柱壳正面、侧面和背面的压力场，得到了冲击波压力时程曲线。

圆柱壳体在水下爆炸冲击波作用下的压力空间分布情况的认识尚不明确。本文基于空气动力学理论，从膨胀波与激波的形成、波在膨胀波与激波中流动过程出发，来解释冲击波绕射的形成，并利用ALE流固耦合算法，探究圆柱壳结构周围压力分布特点与绕射波特性。

1 数值分析模型

1.1 物理模型

本文所采用的物理模型以圆柱壳中心为坐标原点建立模型，圆柱壳半径为R、厚度为h的薄圆柱壳，壳体内部为空气，壳体外部为水，水流场中有TNT当量为W的药包，药包中心与壳体间的距离为r，如图1所示。

图1 物理模型Fig.1 Physical model

本文采用LS-DYNA中的任意拉格朗日欧拉(ALE)算法实现冲击波与圆柱壳的相互作用模拟。依据上述物理模型并结合文献[6]中的实验模型尺寸，建立长1.35 m，直径1 m，板厚为10 mm的圆柱壳体。模型横剖面图如图2所示，壳体内部为空气，外部为水，流场网格大小为0.04 m，结构网格大小为0.02 m。流场边界采用无反射边界条件。基于计算模型的对称性，在壳体上设置A、B、C、D、E5个测点。

1.2 模型参数设置

1)结构材料特性。

本文所采用高强度低碳钢，考虑材料应变率敏感性的本构方程有很多，本次圆柱壳体采用与实验数据吻合较好的Cowper-Symonds模型来描述材料的应变率效应，动态屈服强度：

(1)

式中：σd为动态屈服强度；σ0为屈服极限；E为杨氏模量；ET为切线硬化模量；εP为有效塑性应变；εd为等效塑性应变；D、δ为应变率参数，可由试验获得，文献[7-10]均提到了低碳钢有实验得到的与应变率相关的参数：D=40.4 s-1，δ=5。材料相关参数见表1。其中，ρ1为钢材密度；μ为泊松比。

图2 计算模型Fig.2 Computation model

表1 模型参数Table 1 Model parameters

2)空气与水状态方程。

本次采用流体常用模型MAT _null模拟流场本构关系，空气状态方程采用EOS_linear_polynormal方程，在LS-DYNA程序中，水的状态方程采用Grüneisen状态方程[11]，相关参数[12]见表2、表3。

表2 空气状态方程参数Table 2 Air flow field state equation parameters

表中ρ2代表空气密度，E0代表初始内能，V0代表初始相对体积，C0～C6为方程系数，线性多项式状态方程：

P=C0+C1ζ+C2ζ2+C3ζ3+E′(C4+C5ζ+C6ζ2)

(2)

式中：ζ为当前密度与初始密度之比；E′为当前内能。

表3 水状态方程参数Table 3 Water flow field state equation parameters

线性多项式状态方程见式(3):

(γ0+αμ)E1

(3)

式中：P为压力；μ=ρ/ρ0-1；ρ3为水密度；E1为单位体积内能；S1、S2、S3、γ0、α为常数。水的密度约为1 000 kg/m3，波在水中的传播速度c约为1 500 m/s，水的运动粘性系数约为ν=1.13×10-3N·s/m2。

3)炸药状态方程。

H. Hamashima等[13]通过水下爆炸实验拟合了JWL状态方程，使其能够更好地反应炸药膨胀过程中的二次反应和约束条件。对于TNT炸药，采用JWL[14]状态方程：

(4)

式中：V当前相对体积，A1、A2、r1、r2为常数，ω为比热，各参数的取值表4。

表4 JWL状态方程参数Table 4 JWL state equation parameters

1.3 压力对比分析

1.3.1 自由场压力对比分析

为验证爆炸冲击波载荷的准确性，对自由场中的水下爆炸冲击波进行数值模拟。计算中TNT当量为1 kg，药包初始半径R0=0.053 m，采用中心点爆的形式。压力测点与爆心的距离进行无量纲化r→r/R0，测点距爆心的距离r=10，r=13，r=16，r=19，计算域为边长为2.5 m的正方体区域，由于流场对称，所以只取1/2流场进行计算。流场边界采用无反射边界条件，流场剖面采用对称边界条件，流场网格大小为0.04 m。

图3中给出了自由场中水下爆炸计算的数值模型和典型时刻自由场中的压力云图。图4中给出了r=10，r=13，r=16，r=19处测点的冲击波压力与Zamyshlyayev[15]公式值的对比曲线，从图中可得数值模拟结果衰减过程与经验公式吻合良好，压力曲线在峰值过后出现多次压力震荡，Cole[16]给出的解释是由于装药中不同位置炸药颗粒(对应于数值计算模型中不同炸药单元)爆炸后冲击波叠加引起的[17]。

图3 计算模型与典型时刻压力云图Fig.3 Computation model and pressure cloud chart

1.3.2 实验值与数值计算值的比较

根据文献[6]中的实验，设置1 kg药量，4 m爆距的工况下进行数值模拟计算，得到圆柱壳迎爆面中心、上表面中心、背爆面中心的压力峰值，与实验值比较，如表5所示。从表中可以看出，实验值与仿真值的误差在10%左右。

表5 压力峰值对比表Table 5 Peak pressure contrast

2 冲击波绕射特性分析

2.1 绕射波形成原理

关于传统绕射波现象的研究虽然比较普遍，且主要集中在空气场中。有关水下爆炸冲击波作用下结构毁伤模式和特征的研究虽然已有不少，但水下强冲击波与圆柱体耦合绕射波形成原理分析较少，也没有形成相关理论体系。当水下爆炸冲击波在传播过程中遇到壁面时，其传播方向会发生改变，但传播方向是如何改变的，其中的作用原理是什么？本文基于空气动力学理论，从膨胀波与激波的形成、波在膨胀波与激波中流动过程出发，来解释冲击波绕射的形成。

超声速气流沿外折微小角度的壁面流动时，产生弱膨胀波，气流经过此波后速度增大，压强、密度、温度降低。二维无粘性超声速气流绕折角为dθ的外凹角流动，如图5所示。在这点O之前，气流平行于壁面AO，流至点O，壁面向外转折dθ，流动空间扩大，气流膨胀，并伴随着压强、温度减小。压强变化的扰动从点O向外传播，形成一道弱膨胀波，它和波前气流方向之间的夹角μ1只取决于马赫数Ma1，即

(5)

图4 爆炸冲击波数值压力与经验公式值对比Fig.4 Comparisons numerical value and empirical formulas

图5 微小外折角产生的膨胀波Fig.5 Expansion wave produced by small external bend angle

二维无粘性超声速绕外凹曲面的流动如图6所示。曲面的弯曲部分O1On可看成由无限多个连续转折无限小折角构成。从弯曲部分上的每一点都发出一条弱膨胀波，这无限多条弱膨胀波构成“连续膨胀波”，气流经过这些膨胀波后，马赫数由Ma1增大到Ma2，方向由平行于壁面AO1折转到平行于壁面OnB，这一过程是渐变等熵过程。

根据极限的概念，一个连续的圆柱形拱壳可以看成是由无数段外折的微元折壁所组成的，如图7所示。当爆炸冲击波绕壁面传播时，圆柱壳壁面上的每一个点都相当于一个折点，因此，每一点都将发出一条弱膨胀波，冲击波每经过一道膨胀波，参数值有一个微小的变化，通过一段膨胀波后，参数值及折角发生一个有限量的变化，从而形成了传播方向发生改变的绕射波。

图6 超声速绕外凸曲面的流动Fig.6 Supersonic flow around convex surfaces

图7 绕射波的形成Fig.7 Formation of diffracted waves

2.2 爆炸冲击波的传播过程

采用ALE流固耦合法对模型进行数值仿真计算。计算得到的流场压力云图如图8所示。装药起爆后，当冲击波作用在壳体后，由于壳体内部介质为空气，大部分能量以冲击波的形式返回到流场中，形成反射冲击波，如图8(a)所示；反射冲击波向爆心方向传播，当传播至气泡表面后，由于气泡内部介质的阻抗低于气泡外部介质的阻抗，又有一部分能量以稀疏波的形式反射回到流场中，另一部分则透射到气泡内部，气泡反射回的稀疏波与壳体反射的稀疏波相遇后会相互叠加，使压力迅速降低，如图(b)；初始冲击波继续向前传播，当冲击波绕过圆柱壳顶部和底部时，形成膨胀波，冲击波的传播方向发生改变，绕过临界位置向背爆面传播，如图8(c)、(d)所示。

图8 流场压力云图Fig.8 Pressure distribution map of flow field

2.3 近壁面压力分析

2.3.1 近壁面压力空间分布

取图2中所示的5个近壁面测点为研究对象。工况1～4的壳体厚度为10 mm，对爆距r进行无量纲化处理，r→r/R(R为壳体半径)，爆距分别为0.5、1、2、3；工况5～8在的爆距为2，壳体厚度分别为5、10、15、20 mm。

图9、10为圆柱壳表面测点A～E的压力在不同工况下的时程曲线，测点位置如图2所示。当初始爆炸冲击波传播到A点后，A点出现第1次压力峰值，随后压力迅速衰减，由于壳体中应力波的反射形成了第2个压力峰值；B点处压力峰值降低，压力变化规律与A点相似；测点C位于临界点，压力峰值降低，脉宽增大；D、E2点压力变化规律与迎爆面测点压力变化规律差异较大，绕射冲击波的到达形成的压力值较低，脉宽增加。通过图9和图10的对比可以发现，板厚主要影响应力波的透射规律，各点的载荷时历特征受爆距参数影响更加明显。

图9 不同爆距下各测点压力时程曲线Fig.9 Time-history curve of pressure at different explosion distances

图10 不同板厚下各测点压力时程曲线Fig.10 Time-history curve of pressure at different thicknesses

2.3.2 近壁面压力峰值特征

为探究冲击波近壁面压力的峰值在圆柱壳周围的分布特征，定义无量纲相对压力峰值λPi=Pmi/PmA，其中Pmi表示第i位置处的压力峰值。λP变化趋势如图11所示。

图11 λP值折线图Fig.11 The broken-line graph of λP values

从图11中可以看出，λPB的值最大，临界点和背爆面测点λPC、λPD、λPE的值依次减小；临界点C处λPC与迎爆面处的λPB值相差较大，与背爆面处的λPD、λPE值相差较小；随着爆距与板厚的增加，λPB-λPE值也随之增大。不管爆距与结构形式如何改变，背爆面λPD、λPE的值始终非常接近。

2.3.3 脉宽变化特征

压力峰值所代表的是瞬时的压力，而脉宽代表了冲击波作用的时间，同样是一个重要的参量，定义脉宽是指冲击波从峰值衰减到峰值的1/e所需的时间，定义无量纲相对脉宽λτi=τi/τA，其中τi表示第i位置处的脉宽。图12给出了各工况下每个测点的λτi=τi/τA。背爆面测点处的冲击波脉宽明显比迎爆面测点处的脉宽大。无论爆距与板厚如何改变，背爆面两测点的λτD-E始终非常接近；并且随着爆距与板厚的增加，λτD-E与界点处λτC的值趋于一致。迎爆测点λτB和临界点λτC几乎不发生变化。