水下爆炸载荷作用下变截面结构中应力波的传播特性研究

孟子飞，徐思博，刘文韬，陈海龙

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

0 引 言

近年来，国内外众多学者对水下爆炸结构毁伤的力学机理进行了系统的研究[1-2]，但对于在水下爆炸载荷作用下应力波在舰船结构中的传播特性研究甚少。目前，对于应力波在变截面结构中传递规律的研究主要以变截面杆为主，其中，王礼立等[3]运用一维杆中应力波的初等理论对应力波在锥形放大器中传播特性进行了理论分析，指出了Leftheris 计算公式的错误之处并与实验值进行对比，验证了理论值的正确性，但计算得到的峰值应力比实测值偏大。周光泉、刘孝敏等[4-5]在此基础上运用二维理论研究该问题并取得了较好的结果。但是，数值模拟与理论推导相比有着特殊的优势，如果模型比较复杂并且载荷呈现强非线性时则难以利用理论推导出某处的应力状态。另外，近年来随着SHPB 实验装置研究领域的拓宽，其传统的2 个基本假定（一维应力假定和均匀性假定）均受到挑战[6-8]，而数值模拟对于优化实验设计有着明显的优势，一方面数值模拟可以为实验提供预测，如果数值模拟精度足够的话完全可以代替实验；另一方面在数值模拟中可以对应力波在不同材料、不同截面压杆中的传播特性进行系统的数值分析[9]。此外，在水下爆炸冲击载荷作用下研究变截面杆以及舷侧外板结构中应力波传播的问题仍比较少见。因此，本文尝试结合应力波理论，同时借助有限元的CEL 方法来探讨在水下爆炸载荷作用下弹性波在变截面杆以及舷侧外板结构中的传播特性问题，并且利用冲击响应谱分析方法探讨舷侧外板结构中的响应规律，为舰船遭遇水下爆炸后结构中应力波的传播特性以及冲击响应规律研究提供参考。

1 变截面杆中的应力波理论

本节在材料为线弹性的基础上利用冲量定理、牛顿第三定律分析了弹性波在变截面杆中的反射和透射规律。

在如图1 所示坐标系中研究一等截面均匀杆微元段在应力波下的纵向运动，其中A0为截面面积，为密度，F 为力。

图 1 等截面均匀杆微元段示意图Fig. 1 The sketch of a uniform cross-section bar

在推导控制方程之前，先作2 个基本假设：

1）杆在变形时横截面保持为平面；

2）应力和应变之间遵循胡克定律。

规定应力和应变均以拉为正，质点速度以X 轴正向为正。

根据冲量定理，有

进而

将式（2）代入式（3）可得

联立式（4）和式（5）得

2 水下爆炸流固耦合作用的CEL 方法

2.1 CEL 方法介绍

基于网格的数值方法有拉格朗日方法和欧拉方法，这2 种方法各有优缺点，拉格朗日方法能够很好描述物质结构的运动，而欧拉方法能够较好处理流体大变形的问题。由于拉格朗日方法和欧拉方法具有互补的性质，于是发展形成了拉格朗日方法和欧拉方法相互耦合的CEL 方法（耦合的欧拉-拉格朗日方法）。在CEL 算法中，欧拉网格和拉格朗日网格相互独立。欧拉网格固定不动，材料在网格里流动，材料运动的同时它的各个物理量被映射到欧拉网格中，欧拉网格是独立于空间位置的“监测点”实时追踪流体的通量。而在拉格朗日网格中，网格节点固定在材料内部，网格随着材料的变形而变形，质量、动量和能量随着网格单元的运动而传递。在CEL计算时，结构利用拉格朗日方法求解，流体采用欧拉方法求解。在流体和结构相互耦合的界面上，拉格朗日单元受到欧拉材料给予的压力载荷，欧拉单元接受拉格朗日单元施加的几何约束，二者计算信息交换通过2 组网格之间的映射实现。基于以上原理，CEL 算法在解决物体的大位移，如碰撞、流体动力学及流体与固体之间的相互作用时，具有明显优势[11]。

2.2 控制方程

在CEL 方法中，拉格朗日数值方法的控制方程[12]如下：

根据质量守恒有

根据动量守恒有

根据能量守恒有

式中： ρ，v，σ，b，e分别为拉格朗日流体介质的密度，速度，应力，质量力和内能。

在CEL 方法中，欧拉数值方法的控制方程[12]如下：

欧拉算法中，材料在固定的网格里流动，根据质量守恒定律，有连续性方程

根据动量守恒有

根据能量守恒有

2.3 欧拉体积分数法

在CEL 中，材料在每一个增量步都会计算每个单元材料的分布以此来描述流体的变形状态，而材料的初始分布情况在Abaqus 中利用“体积分数工具”进行计算。该工具对欧拉体与参考体之间进行布尔运算来确定材料的分布。图2 中利用“体积分数工具”将T 型材作为参考体映射到欧拉体后计算得到的欧拉体积分数示意图，图中灰色部分代表在欧拉网格中填充“水”材料，数字代表“水”在欧拉单元中所占的体积分数，白色部分代表T 型材结构其材料的敷设在拉格朗日体中进行。后处理中通过EVF 选项卡来观察流体的运动情况。

图 2 欧拉体积分数示意图Fig. 2 The sketch of eulerian volume fraction

2.4 罚函数接触方法

在CEL 方法中，欧拉域与拉格朗日域之间的接触是通过罚函数接触方法实现的。拉格朗日单元边和面上建立节点，欧拉材料表面建立锚点，二者通过“硬”接触-压力过盈行为建立耦合关系[13]。这种方法允许欧拉材料和拉格朗日体之间存在一些微小的渗透。节点与锚点之间的接触力Fp正比于二者的渗透距离d[14]。

式中：k 为罚函数刚度系数，取决于拉格朗日和欧拉材料的属性。

2.5 CEL 方法验证

本节采用自由场水下爆炸模型，以仿真值与经验公式值对比的方式对CEL 方法进行验证，具体的几何模型如图3 所示。

球形T N T 的药量为0.0 5 5 g，装药密度1 600 kg/m3，爆速 。爆轰产物的压力与爆轰产物相对比容和比内能之间的关系釆用JWL 状态方程表示[13]。方程中参数具体值见表1。

水的密度ρ=1 000 kg/m3，粘性系数为0.001，采用US-UP形式的EOS 状态方程，其中，水的声速C0=1 480 m/s。

图 3 自由场水下爆炸几何模型示意图Fig. 3 The geometry model for underwater explosion of free field

表 1 JWL 方程的参数Tab. 1 The parameter of JWL equation

为减小边界对冲击波影响，计算欧拉域的所有边界都设定为流出无反射边界，在无边界效应影响的水域进行爆炸，其产生冲击波的波阵面压力峰值pm可由Cole 和Zamyshlyayev 给出的经验公式[15-16]表示。图4为在不同距离处CEL 计算压力峰值与经验公式计算压力峰值的对比。从图中可以发现，CEL 方法计算出来的结果与Cole 和Zamyshlyayev 给出的经验公式计算结果相似，准确度较高。因此，采用CEL 方法可以充分模拟水下爆炸冲击载荷。

图 4 不同距离处的压力峰值对比Fig. 4 Peak pressure contrast at different distances

3 水下爆炸载荷作用下变截面杆中应力波的传播特性

下面采用CEL 方法模拟出水下爆炸冲击载荷，然后通过有限元方法模拟冲击波在变截面杆中的传播问题。此问题涉及液固相互作用并且水下爆炸载荷呈现非线性难以利用理论公式得到解析解，本节运用CEL方法来处理流固耦合问题以期达到研究变截面杆在受到水下爆炸载荷下弹性波传播特性的目的。

本节建立1 根两段等长L 的变截面杆，L 为200 mm，杆左端直径D 为4 mm，杆的截面比为2∶1，杆的材料为普通钢（不考虑塑性），杆的两端均为自由端，在距杆左端50 mm 处建立和第2.5 节一样的TNT 炸药模型，最后建立一足够大的Euler 域以保证TNT 炸药具有足够的传播空间，Euler 域四周均设置无反射边界条件，具体模型示意图如图5 和图6 所示。

图 5 变截面杆几何示意图Fig. 5 The geometry model of a variable section bar

图7～图9 给出了爆炸冲击波从开始进入到杆的左端、传递到变截面附近、传递到杆右端附近3 个典型位置图。当时间s 时，提取出入射波作用下杆左端到右端的应力随距离传递曲线和速度随距离传递曲线，计算距离从0～0.2 m 对应的应力峰值作为速度峰值作为 vI=0.28 m/s。当时间s 时，提取出透射波作用下杆左端到右端的应力随距离传递曲线和速度随距离传递曲线，计算距离从0.2～0.4 m 对应的应力峰值作为，速度峰值作为m/s。各曲线如图10 所示，与理论解进行对比如表2 所示。可以发现变截面杆在受到水下爆炸载荷作用下弹性波规律十分接近理论解，并且可以看出当应力波由大截面向小截面传播时，透射波扰动强度大于入射波扰动，这种结构起到“聚波”作用。由于在t=1.32×10-4s 时，杆左端到杆变截面处存在反射波与部分入射波的叠加，所以难以提取反射波作用下杆的应力传递曲线和速度传递曲线，在此不予处理。

图 7 冲击波进入杆左端典型位置图Fig. 7 The shock wave enters the left of the bar

图 8 冲击波传递到变截面附近位置图Fig. 8 The shock wave enters the center of the bar

图 9 冲击波传递到杆右端附近位置图Fig. 9 The shock wave enters the right of the bar

4 水下爆炸载荷作用下舷侧外板中应力波的传播特性

4.1 模型建立

本节借助CEL 方法以期达到研究舷侧外板在受到水下爆炸载荷下弹性波传播特性的目的。此问题与上节变截面杆模型相比，不仅涉及液固相互作用而且模型更为复杂，无法利用理论公式得到解析解。

图 10 杆的应力和速度随距离变化曲线Fig. 10 Stress and velocity change with distance

表 2 仿真值与理论公式值对比Tab. 2 Comparison of simulation and theoretical values

图 11 舷侧外板遭遇水下爆炸仿真模型Fig. 11 The simulation model of side plating underunderwater explosion

4.2 应力波从舷侧外板到T 型材腹板传播的特性及冲击响应谱分析

采用频谱分析方法和冲击响应谱分析方法，探究应力波从舷侧外板结构往T 型材腹板中的传播特性。分别提取出#1，#2，#7，#8 四个测点的加速度随时间变化的数据，并进行快速傅里叶变换得到对应的频谱曲线，如图12 所示。从图中可以看出，在[500,2 000] Hz 较高频率范围内的应力波从舷侧外板经过第1 个T 型材后传递到T 型材腹板处的强度大为减小，而[0, 500] Hz 内的应力波强度则相差不大。这说明应力波在经过第1 个从小到大截面时滤掉了大部分高频信号。

图 12 各测点频段分析Fig. 12 Analysis on frequency band at various gauging points

将以上测点的加速度时历数据进行转换得到冲击响应谱曲线，如图13 所示。可以发现，对于[500,2 000] Hz 较高频率范围内的冲击响应经过第1 个T 型材后呈大幅衰减趋势。这说明舷侧外板在受到水下爆炸载荷作用下[500, 2 000] Hz 内的T 型材腹板的结构响应与外板结构响应相比大为减小。而[0, 500] Hz 内的响应相对来说变化则较为平缓。总体来看，从舷侧外板到第1 个T 型材，低频段的结构响应谱值两者相差不大，而在高频段外板的谱值远高于T 型材腹板中的谱值。

图 13 各测点处冲击响应谱曲线Fig. 13 The curve of shock response spectrum at various gauging points

分别提取出#9，#10 测点的加速度随时间变化的数据，并进行快速傅里叶变换，得到对应的频谱曲线，并与#3，#4 进行对比，如图14 所示。可以看出，在[0, 500] Hz 的应力波从舷侧外板经过第2 个T 型材后传递到T 型材腹板处的强度大为减小，而在其他频段中腹板的应力波强度与外板相比相差不大。这说明当应力波再次经过从小到大截面形式的变化时滤掉了大部分低频信号。

图 14 各测点频段分析Fig. 14 Analysis on frequency band at various gauging points

将以上测点的加速度时历数据进行转换得到冲击响应谱曲线，如图15 所示。可以看出，冲击响应规律与经过第1 个T 型材后的规律不同，舷侧外板在受到水下爆炸载荷作用下，只有[1 300, 1 800] Hz 内的结构响应才大于腹板的响应，在其他频段两者则呈波动交替趋势，幅值相差不大。

图 15 各测点处冲击响应谱曲线Fig. 15 Shock response spectrum at various gauging points

结合上述描述可以得出，当应力波再次经过从小到大截面形式的变化时，舷侧外板产生的中高频段冲击响应大于T 型材腹板中的响应，而其他频段的响应相差不大。

4.3 应力波在舷侧外板中的传播特性及冲击响应谱分析

本节采用上节相同的分析方法，进一步探究应力波在舷侧外板结构中的传播特性。

分别提取出#1～#6 六个测点的加速度随时间变化的数据，并进行快速傅里叶变换，得到对应的频谱曲线，如图16 所示。

通过图16（a）与图16（b）两组测点相互对比可以发现，在[500, 2 000] Hz 内#3，#4 测点的幅值比#1，#2 测点小，这说明在该频率范围内的应力波经过第1 个T 型材后强度呈大幅衰减趋势，衰减幅度约60%，而#3，#4 测点在[0, 500] Hz 内的幅值比对应的#1，#2 测点大，这说明在该频率范围内的应力波经过第1 个T 型材后强度略微增强。结合这两点，说明当应力波经过第1 个T 型材后滤掉了大部分高频信号，而增强了低频信号，前者是由于应力波在经过第1 个从小截面到大截面形式的变化时会滤掉高频信号的特性所导致，后者是因为应力波从大截面进入小截面其强度会增强的原理所导致。

图 16 各测点频段分析Fig. 16 Analysis on frequency band at various gauging points

从图16（c）与图16（d）可以发现，#5，#6 测点在[0, 500] Hz 内的幅值比对应的#3，#4 测点小，这说明在该频率范围内的应力波经过第2 个T 型材后强度呈衰减趋势，衰减幅度约50%，而#5，#6 测点在[500,2 000] Hz 内的幅值相对于#3，#4 测点来说衰减幅度并不大。基于此，这说明应力波再次经过相同截面形式的T 型材后会滤掉大部分低频信号，这与4.2 节的结论相吻合。

将以上测点的加速度时历数据进行转换得到其对应的冲击响应谱曲线，如图17 所示。可以看出，#1，#3，#5 测点处冲击响应谱曲线与#2，#4，#6 测点处冲击响应谱曲线规律相似。为方便起见，下面以#2，#4，#6 测点处冲击响应谱曲线规律进行阐述。在[0,2 000] Hz 内，#2，#4，#6 测点的谱峰值依次减小，说明舷侧外板在受到水下爆炸载荷作用下，在[0,2 000] Hz 内的结构响应峰值在经过2 个T 型材后呈衰减趋势。具体说来，在[0, 500] Hz 内#2 测点的谱峰值小于#4，则说明舷侧外板在受到水下爆炸载荷作用下[0, 500] Hz 内的结构响应在经过第1 个T 型材后呈放大趋势，而#6 测点的谱峰值小于#4 测点，这说明舷侧外板在受到水下爆炸载荷作用下[0, 500] Hz 内的结构响应在经过第2 个T 型材后呈衰减趋势。对于[500,2 000] Hz 较高频率范围内的冲击响应经过2 个T 型材后整体呈衰减趋势。

综上所述，舷侧外板中低频段的冲击响应谱峰值沿着从下往上的方向先升高再下降而高频段的谱峰值则依次下降。

图 17 各测点处冲击响应谱曲线Fig. 17 The curve of shock response spectrum at various gauging points

5 结 语

本文基于应力波理论，结合CEL 方法，首先研究了在水下爆炸载荷作用下变截面杆结构中应力波的传播特性，与理论公式对比验证了数值精度，然后进一步探讨了舷侧外板结构在遭遇水下爆炸载荷后的应力波传播特性，得到如下结论：

1）当应力波由大截面向小截面传播时，透射波扰动强度大于入射波扰动，这种结构起到“聚波”作用。

2）应力波在经过第1 个从小到大截面形式的变化时滤掉了大部分高频信号，当应力波再次经过从小到大截面形式的变化时会进一步滤掉响应信号。

3）在舷侧外板到T 型材腹板这种从小到大截面形式变化的结构中，激励产生的低频段结构响应谱值，两者之间相差不大，而外板的高频段谱值远高于T 型材腹板中的谱值。当应力波再次经过这种截面形式的变化时，舷侧外板产生的中高频段冲击响应大于T 型材腹板中的响应，而其他频段的响应相差不大。

4）在含有2 个T 型材的舷侧外板中，低频段的冲击响应谱峰值沿着从下往上的方向先升高再下降而高频段的谱峰值则依次下降。

以上结论对研究船舶与海洋工程领域的型材、纵骨、纵桁等变截面结构中应力波传播特性以及后续研究复杂舰船结构的冲击响应规律具有一定的参考意义。