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

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

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

1 数值分析模型 1.1 物理模型

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

本文采用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模型来描述材料的应变率效应，动态屈服强度：

$\sigma_{d}=\left(\sigma_{0}+\frac{E E_{T} \varepsilon_{P}}{E-E_{T}}\right)\left[1+\left(\frac{\varepsilon_{d}}{D}\right)^{\frac{1}{\delta}}\right] $

(1)

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

2) 空气与水状态方程。

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

表中ρ₂代表空气密度，E₀代表初始内能，V₀代表初始相对体积，C₀~C₆为方程系数，线性多项式状态方程：

$P=C_{0}+C_{1} \zeta+C_{2} \zeta^{2}+C_{3} \zeta^{3}+E^{\prime}\left(C_{4}+C_{5} \zeta+C_{6} \zeta^{2}\right) $

(2)

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

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

$\begin{array}{*{20}{c}} {P = \frac{{{\rho _3}{c^2}\mu \left[ {1 + \left( {1 - \frac{{{\gamma _0}}}{2}} \right)\mu - \frac{\alpha }{2}{\mu ^2}} \right]}}{{{{\left[ {1 - \left( {{S_1} - 1} \right)\mu - {S_2}\frac{{{\mu ^2}}}{{1 + \mu }} - {S_3}\frac{{{\mu ^3}}}{{{{(1 + \mu )}^2}}}} \right]}^2}}} + }\\ {\left( {{\gamma _0} + \alpha \mu } \right){E_1}} \end{array} $

(3)

式中：P为压力；μ=ρ/ρ₀-1；ρ₃为水密度；E₁为单位体积内能；S₁、S₂、S₃、γ₀、α为常数。水的密度约为1 000 kg/m³，波在水中的传播速度c约为1 500 m/s，水的运动粘性系数约为ν=1.13×10^-3 N·s/m²。

3) 炸药状态方程。

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

$p=A_{1}\left(1-\frac{\omega}{r_{1} V}\right) e^{-r_{1} V}+A_{2}\left(1-\frac{\omega}{r_{2} V}\right) \mathrm{e}^{-r_{2} V}+\frac{\omega e}{V} $

(4)

式中：V当前相对体积，A₁、A₂、r₁、r₂为常数，ω为比热，各参数的取值表 4。

1.3 压力对比分析 1.3.1 自由场压力对比分析

为验证爆炸冲击波载荷的准确性，对自由场中的水下爆炸冲击波进行数值模拟。计算中TNT当量为1 kg，药包初始半径R₀=0.053 m，采用中心点爆的形式。压力测点与爆心的距离进行无量纲化r→r/R₀，测点距爆心的距离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]。

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

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

2 冲击波绕射特性分析 2.1 绕射波形成原理

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

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

$\mu_{1}=\arcsin \frac{1}{M a_{1}} $

(5)

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

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

2.2 爆炸冲击波的传播过程

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

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、E 2点压力变化规律与迎爆面测点压力变化规律差异较大，绕射冲击波的到达形成的压力值较低，脉宽增加。通过图 9和图 10的对比可以发现，板厚主要影响应力波的透射规律，各点的载荷时历特征受爆距参数影响更加明显。

2.3.2 近壁面压力峰值特征

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

从图 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几乎不发生变化。

2.3.4 背爆面近壁面相邻位置压力分布特征

经上述分析发现背爆面与迎爆面处的压力分布特征存在较大差异，本节分析背爆面相邻测点处的压力分布特征。将D、E 2处的峰值之比和脉宽之比在不同爆距与板厚下的变化趋势绘于图 13中。

可以发现，随着爆距的增大，P_mD/P_mE从1.16降低至0.99，在爆距为3时，测点E的压力峰值超过了测点D的压力峰值，脉宽之比τ_D/τ_E从0.79降至0.97。

3 结论

1) 水下装药爆炸后释放大量能量，以冲击波的形式向外传播，当冲击波传播到结构物表面时，出现压力峰值；一部分冲击波会反射回流场形成反射冲击波，入射波传播到临界点时会改变其传播方向，形成绕射波到达平板背爆面。

2) 迎爆面测点处压力迅速达到峰值，压力峰值大但脉宽较小；圆柱壳顶部及背爆面处测点压力上升缓慢，压力值小但脉宽较大。迎爆面的压力峰值相比与背爆面大致相差一个数量级，而背爆面的冲击波脉宽时间能达到迎爆面的数倍。

3) 迎爆面测点的压力分布特征与背爆面测点的压力分布特征存在明显差异，临界点的压力特征与背爆面压力特征更为相似。无论爆距与板厚如何改变，背爆面各处的压力峰值与脉宽不随空间位置的变化而变化。

4) 随着爆距和板厚的增加，背爆面的冲击波脉宽λ_τD-E与临界点处的冲击波脉宽λ_τC的值趋于一致。