谢富佩，徐 绯，曾 卓，周中玉，谷卓伟

（1. 西北工业大学航空学院计算力学与工程应用研究所，陕西 西安 710072；2. 中国航空发动机集团有限公司湖南动力机械研究所，湖南 株洲 412002；3. 中国工程物理研究院流体物理研究所，四川 绵阳 621999）

复合圆柱壳体结构具有高比刚度、高比吸能以及轻质等优点，在航空、航海、汽车工业及军事等领域有广泛应用。在强磁场物理、材料合成等领域，使用了一种利用炸药爆轰驱动复合圆柱壳内部预置磁场以得到高能磁场的装置，通常将此装置称为柱面内爆磁通量压缩发生器(explosive magnetic flux compression generator)，简称MC-1 装置[1]。与常见的纤维编织复合圆柱壳不同，由于考虑到复合圆柱壳通电后产生并压缩磁场的要求，此类复合结构中铜线纤维是连续平行排列的；在炸药爆轰产生的冲击载荷作用下，MC-1 装置中复合圆柱壳容易发生界面失稳和扰动增长而导致结构失稳[2]；另外，结构也可能在冲击压缩过程中发生较大的塑性变形而失效，而复合圆柱壳的结构响应及稳定性是影响MC-1 装置性能的重要因素。目前MC-1 装置在高能量密度物理实验技术研究领域应用较为广泛[3]，因此对复合圆柱壳体结构在冲击压缩下的结构响应模拟及稳定性研究有十分重要的意义。

目前国内外对于MC-1 装置中的复合圆柱壳在冲击压缩下的数值模拟研究较少。俄罗斯实验物理科学研究所于2015 年公布了几种MC-1 装置的设计参数，但未公开其参数设计细节；陆禹等[4]使用了引入解析形式状态方程组的一维内爆磁流体力学编码MC11D，对MC-1 装置工作过程进行了一维条件下的数值模拟；张春波等[5]利用AUTODYN 软件二次开发接口建立了内爆压缩多层密绕螺线管复合圆柱壳三维模型，并实现了周期性边界条件，使用光滑粒子流体动力学方法对其动态响应及不稳定性开展了数值模拟；赵士操等[6]使用SPH 方法构建了纤维织物的细观尺度计算模型，使用AUTODYN 软件对超高速碰撞下纤维增强复合材料模型进行了分析计算；刘军等[7]对爆轰驱动内壁刻有正弦扰动的金属钢壳与内部硅橡胶界面产生不稳定性问题进行了数值模拟。由以上研究可以看出，对炸药爆轰驱动作用下的复合圆柱壳动力学响应及稳定性数值模拟问题，由于涉及炸药爆炸、高速冲击等过程，传统网格法面临大变形、高度非均匀性、运动物质交界面及边界处理等诸多问题，而粒子法又有计算耗时较长的问题，因此对复合圆柱壳细节模型的数值模拟常存在较大困难，导致现有的研究缺乏对复合圆柱壳细节模型的分析手段，且实验对结构失稳现象未提出合理的失稳判据，复合圆柱壳稳定性的影响因素难以评估。

基于以上缺陷，本文以MC-1 装置中的复合圆柱壳作为研究对象，使用光滑粒子动力学-有限元方法（smooth particle hydrodynamic-finite element method, SPH-FEM）耦合算法，构建复合圆柱壳二维细节爆轰模型并对其在炸药爆轰冲击下的结构响应进行分析计算，提出一种判断圆柱壳结构在冲击压缩下的稳定性评价方式，分析影响MC-1 装置性能的因素，以期为MC-1 装置设计及相关实验研究提供参考。

1 数值计算模型

1.1 模型介绍

本文计算模型的结构参数参照俄罗斯实验物理科学研究所核中心于2015 年公布的实验参数，装置结构如图1(a)所示，复合圆柱壳外侧包裹着一层环形炸药，炸药外侧均匀布置48 个起爆点，在同一时刻起爆。装置中的复合圆柱壳在制备时首先使用一定直径的漆包铜导线按照规定的角度螺旋绕制成多匝多层的螺线管内筒体，将绕制的螺线管表面使用环氧混合剂固定并形成一层绝缘层，然后将从螺线管伸出的铜线覆盖在绝缘层表面，覆盖的铜线方向与复合圆柱壳轴线方向平行，最后在其表面使用环氧混合剂凝固。因此MC-1 装置中复合圆柱壳由内向外分别为铜线层、内环氧层、铜线折返层及外环氧层，如图1(b)所示。图1(c)为公布的实验装置的结构参数示意图，考虑到理想情况下铜线层及铜线折返层的厚度为铜线直径的整数倍，因此计算时将铜线层的厚度调整为2 mm，进行模拟计算时装置的参数如图1(d)所示。

图1 复合圆柱壳示意图Fig. 1 Diagram of composite cylindrical shell

1.2 数值模型的建立

对在爆轰驱动作用下复合圆柱壳结构响应的模拟计算问题，由于复合圆柱壳三维模型较为复杂，且计算量较大，而复合圆柱壳结构沿轴线方向的截面形状及尺寸无明显变化，因此本文选用对复合圆柱壳具有代表性的二维截面模型进行计算。考虑到只使用有限元算法时装置内圆柱壳结构可能会出现的网格畸变问题以及只使用SPH 算法可能导致的计算规模较大的问题，本文使用SPH-FEM 耦合算法[8]对爆轰冲击下的复合圆柱壳进行数值模拟，即对装置炸药部分使用FEM 方法划分网格，对本问题的主要研究对象复合圆柱壳部分使用SPH 方法建立模型，由于两种方法都具有明显的拉格朗日特性，因此只需要在软件中设置接触域即可实现数据传递，完成两种算法的耦合。在建立模型时，由于复合圆柱壳结构较为复杂，直接在AUTODYN 软件建模较为困难，本文使用导入外部文件[9]的方式建立模型。首先在CATIA 软件中建立复合圆柱壳细节模型并导入到有限元软件HyperMesh 中划分三角形网格，在划分网格时，网格尺寸即为后续计算的粒子间距，本文将粒子间距取为0.1，导出K 文件后编写C 程序将模型节点格式转化为AUTODYN 软件可识别的SPH 粒子格式，将其存储为Formatted 文件后导入AUTODYN软件，进一步添加有限元炸药，完成模型边界和起爆点等的设置后即可完成建模。图2 为最终在AUTODYN 中建立的15°铜线螺旋角下的复合圆柱壳数值模型。

图2 15°螺旋角下的爆轰驱动复合圆柱壳模型Fig. 2 Detonation driven composite cylindrical shell model with a spiral angle of 15°

在数值模拟中，复合圆柱壳外侧用于提供爆轰压力的炸药部件使用JWL 状态方程[10]描述：

表1 Comp B 炸药主要材料参数Table 1 Main material parameters of Comp B explosive

复合圆柱壳中铜线类型为高导无氧铜，其材料参数[11]如表2 所示，塑性部分使用Johnson-Cook 材料模型，材料的屈服应力可定义为：

表2 铜线主要材料参数Table 2 Main material parameters of copper wire

圆柱壳的环氧部分的材料模型则使用Shock 状态方程描述:

式中：U为激波速度；up为粒子速度；c0和s为常数，参数取值参照AUTODYN 材料库，c0取为2.73 km/s，s取为1.493。

数值计算模型中圆柱壳与炸药之间的接触定义为Lagrange/Lagrange 接触，炸药外侧设置48 个均布起爆点，均在0 µs 时刻同时起爆。

2 结构响应与稳定性分析

2.1 数值模型结果验证与分析

为了验证建立的数值模型计算结构响应的准确性，在AUTODYN 中建立炸药爆轰驱动金属环氧复合套筒模型进行数值模拟，套筒的内径为94 mm，外径为100 mm，沿径向由内向外分别为厚度为1.5 mm钢套筒和1.5 mm 的环氧套筒，外部环形炸药的厚度为55 mm。将计算的套筒内侧半径变化曲线与使用已经过验证的Triangels-MHD 磁流体力学程序[12]计算的结果进行对比，如图3 所示。

图3 套筒内侧半径变化曲线对比Fig. 3 Comparison of radius change curves within the sleeve

根据半径变化对比曲线可以很明显地看出，除了曲线后半段外，使用两种方法计算得到的套筒内壁位移响应曲线大部分阶段都比较接近，而曲线后半段出现一定差异是因为TriAngels-MHD程序主要是基于磁流体力学方程建立的，在求解过程中考虑到了磁力的影响，在压缩后期磁力会由于半径减小而突增，使用AUTODYN 软件求解时则暂未考虑到后期磁力影响。但是由于磁力影响只占套筒压缩过程中的极小一部分，因此可以认为本文使用的二维数值计算模型能较好地模拟圆柱壳在炸药冲击压缩下的结构响应。

而由于金属环氧套筒的计算模型中未引入初始扰动，整个爆轰过程中各个起爆点的起爆时间相同，起爆点分布也较为密集，使得套筒外侧的爆轰载荷沿环向分布较为均匀，因此在压缩过程中套筒无较大的环向扰动，基本保持为一个圆环，如图4(a)所示。同样地，复合圆柱壳的变形也有此特点，以15°铜线螺旋角的复合圆柱壳模型为例，在爆轰作用下其不同时刻的变形如图4(b)所示，从图中可以看出，在压缩过程中，复合圆柱壳的环氧层及铜线折返层有较大的环向变形，而决定装置整体性能的铜线层则在压缩前中期基本保持为环向相对光滑的状态，直至压缩后期才有轻微的环向扰动。

图4 圆柱壳在不同时刻的变形Fig. 4 Deformation patterns of cylindrical shell at different times

2.2 稳定性判据

在MC-1 装置的实际应用中，复合圆柱壳受到爆轰载荷后向内收缩是一种高压高应变率条件下的大变形行为，在此阶段中极易发生界面不稳定增长，而复合圆柱壳失稳是影响MC-1 装置性能的重要因素，因此在数值模拟中需引入一种失稳判据，以对复合圆柱壳在冲击压缩过程中的稳定性进行评价。

首先以金属环氧套筒为研究对象，图5(a)为计算的金属环氧套筒在冲击压缩过程中内壁粒子的速度变化曲线，曲线可大致分为3 段：(1) 炸药爆轰波尚未到达套筒表面，结构尚无响应，速度为0；(2) 爆轰波到达套筒表面，套筒在爆轰波作用下开始向内收缩，内壁粒子速度急剧增加后逐渐趋于平稳；(3) 结构出现扰动增长，内壁粒子速度进一步增加并出现不稳定增长。根据套筒内壁粒子速度变化历史曲线可大致地判断出结构在第3 阶段出现失稳。为了获得更加精确的失稳时间和失稳时对应的压缩率等数据，对第3 阶段的速度曲线进行拟合并对拟合曲线进行求导，图5(a)左上图为对速度曲线第3 阶段拟合曲线求一阶导的导数曲线，导数曲线出现明显转折的时间即为套筒结构的失稳时间。经过曲线处理后，可判断出金属环氧套筒结构的失稳时间大致在15.20 µs 左右，此时套筒的内半径为18.02 mm，整个结构压缩率大致为61.66%。

图5 金属环氧套筒失稳时间判断曲线Fig. 5 Judgment curves of instability time of metal epoxy sleeve

为了证明提出的失稳判据的合理性，沿金属环氧套筒径向由外向内等间距取点，并将其标记为Node 1～Node 5，提取其冲击压缩作用下的速度变化曲线，如图5(b)所示，可以看到在炸药爆轰载荷的作用下，金属环氧套筒径向各层粒子在压缩前期变化规律基本一致，而压缩至15 µs 左右时各层粒子的速度曲线开始出现较为明显的差异，内层粒子的速度增加，而外层粒子速度减小，这种速度差异势必会导致结构径向的失稳，而此失稳时间也与前面提出的使用内壁粒子速度曲线变化判断的失稳时间十分接近，由此可以证明提出的失稳判据较为合理。

将上述失稳判据推广到复合圆柱壳结构，可以发现相同的规律。以15°铜线螺旋角的复合圆柱壳为例，图6(a)为复合圆柱壳结构内壁粒子的速度变化曲线，也有明显的三阶段特征，图6(a)左上图为将复合圆柱壳内壁粒子速度变化曲线的第3 阶段曲线拟合求导后的曲线，可以看到曲线在24.76 µs 时开始有突增的趋势，而此时刻与图6(b)所示的复合圆柱壳沿径向各层粒子速度开始出现明显差异的时间相近，图6(b)中Node 1～Node 11 指沿复合套筒径向由外向内等距标记的11 个粒子。综合以上分析，可以证明本文提出的判断失稳时间的判据对复合圆柱壳结构的稳定性分析有一定的参考价值。在后面的研究中，将继续使用根据结构内壁粒子速度变化历史判断失稳时刻的判据对不同条件下的复合圆柱壳结构的稳定性进行分析讨论。

图6 复合圆柱套筒失稳时间判断曲线Fig. 6 Judgment curves of instability time of composite cylindrical shell

3 影响复合圆柱壳稳定性因素分析

3.1 复合圆柱壳制备缺陷对稳定性的影响

由复合圆柱壳的制备过程可知，铜线折返层的铜线是复合圆柱壳内层的铜线层延伸出来的，因此在理想情况下复合圆柱壳的铜线层和铜线折返层铜线数量应相同。在绕制复合圆柱壳结构的内层铜线层时，铜线层椭圆截面的长轴长b计算公式为

式中：D为铜线直径，θ 为铜线的螺旋角度，同样以15°螺旋角的复合圆柱壳模型为例，内层铜线层铜线截面的长轴长约为0.965 9 mm，则内层每层铜线的数量ni可表示为

式中：Ri为每层铜线的铜线中心距圆柱壳中心的距离，l为铜线层数，按下式计算

式中：h为模型中铜线层的厚度，对本文中计算的铜线直径为0.25 mm 的模型，结合式(5)和式(6)估算出铜线层铜线的数量约为3668 根。

而理想条件下，将模型的铜线折返层全部排满大致需要的铜线数量为5623 根，显然在实际制备复合圆柱壳时并不能将铜线折返层完全由铜线排满，在铜线折返层极有可能存在一定的缺陷。因此建立含有缺陷的复合圆柱壳模型来研究复合圆柱壳在制备过程中的缺陷对其稳定性的影响，具体做法为在复合圆柱壳折返层建模时随机删去若干铜线，本文中含缺陷模型随机在铜线折返层外侧删去100 根铜线，缺陷部分细节如图7(a)所示，模型共有4 处这样的缺陷。根据计算结果可以发现，由于铜线折返层缺陷部分强度较小，在炸药爆轰压缩复合圆柱壳的过程中会导致强度较弱的部分铜线折返层及内环氧层侵入铜线层，并且随着爆轰波的运动，侵入部分的深度增加，影响铜线层的变形形态，进而会影响整个结构的稳定性，不同时刻含缺陷模型的复合圆柱壳局部形态如图7(b)所示，图7(c)为含缺陷模型复合圆柱壳在24.58 µs 的形态，可以看到其结构内壁环向不再保持为圆形，图7(d)为失稳时间判断曲线，含缺陷结构的失稳时间为22.63 µs，装置的稳定性明显降低。以上对含缺陷结构的复合圆柱壳的计算表明：复合圆柱壳结构折返层的制备缺陷会影响结构稳定性。

图7 含缺陷复合圆柱壳变形图及失稳时刻判断曲线Fig. 7 Deformation diagrams of defective composite cylindrical shell and the judgement curve of instability time

3.2 铜线螺旋角对稳定性的影响

复合圆柱壳是非密实、非均匀的结构，铜线的螺旋角度可能是影响其稳定性的重要因素[13]。考虑到在实际应用中复合圆柱壳中铜线螺旋角通常在20°以内，故建立螺旋角度分别为0°和10°的复合圆柱壳模型，不同铜线螺旋角的圆柱壳细节如图8 所示。除螺旋角度外其余参数及条件均与前面的15°螺旋角模型相同，建模时不同角度模型的粒子数也相近。比较不同铜线螺旋角度下复合圆柱壳的结构响应，使用本文第2 节提到的速度判断失稳判据对模型的失稳时间进行判断，结合计算的复合圆柱壳位移变化结果，不同螺旋角模型的失稳时间及失稳时对应的压缩率比较如表3 所示。

表3 不同螺旋角复合圆柱壳失稳时间及压缩率对比Table 3 Comparison of instability time and compression ratio of composite cylindrical shell with different spiral angles

图8 不同螺旋角复合圆柱壳细节模型Fig. 8 Detailed composite cylindrical shell models at different spiral angle

实际上，在建立几何模型时，可以推算出复合圆柱壳铜线部分的截面积S与铜线的螺旋角度无关的结论：

因此，当铜线直径相同时，不同螺旋角度铜线在复合圆柱壳结构中的截面积相同，建立的计算模型中粒子数目相近，因此可得到螺旋角度对复合圆柱壳影响不大的结论。

3.3 铜线直径对稳定性的影响

由式(7) 可知复合圆柱壳铜线部分的截面积与铜线直径有关，另外当铜线直径发生变化时，相应地，复合圆柱壳各层的厚度也有可能发生改变，因此本节讨论铜线直径对复合圆柱壳的影响，主要讨论0.25 和0.50 mm 两种直径，铜线螺旋角度均为15°。当铜线直径为0.50 mm时，可以估算出将2 mm 的铜线层排满需917 根铜线，这些铜线在铜线折返层只需排列一层，图9为0.50 mm 铜线直径下复合圆柱壳模型的结构参数及细节示意图。使用前文提到的失稳判据，可以判断对0.50 mm 铜线直径的复合圆柱壳模型的失稳时间大致为23.62 µs，如图10 所示，此时对应的压缩率为66.41%，比0.25 mm 直径的模型失稳时对应的压缩率小。另外当选用0.50 mm直径的铜线时，由于在复合圆柱壳只有一层铜线折返层，制备时极易出现折返层缺陷而导致装置稳定性变差，图11 为0.50 mm 铜线直径下复合圆柱壳铜线折返层仅有一根铜线缺失时在23.20 µs 时的变形形态，可看到仅缺失一根铜线就会导致复合圆柱壳内壁在压缩过程中有较大的环向扰动。

图9 0.50 mm 铜线直径下复合圆柱壳示意图Fig. 9 Schematic diagrams of composite cylindrical shell under 0.50 mm copper wire diameter

图10 0.50 mm 直径下失稳时间判断曲线Fig. 10 Judgment curve of 0.50 mm model instability time

图11 t=23.20 µs 时含缺陷复合圆柱壳 变形Fig. 11 Deformation of defective cylindrical shell at t=23.20µs

4 结 论

本文利用AUTODYN 软件建立了复合圆柱壳二维细节模型，基于SPH-FEM 耦合算法，对在炸药爆轰冲击压缩下的复合圆柱壳的结构响应过程进行了数值模拟，并根据圆柱壳内壁粒子速度历史提出了一种判断失稳时间的方法。通过对金属环氧套筒爆轰模型的计算，验证了数值模型计算的准确性，另外由复合圆柱壳及金属环氧套筒沿径向不同位置的速度出现明显差异的时间与使用提出的失稳判据判断的失稳时间一致，证明了提出的失稳判据的合理性。

通过对不同数值模型的计算对比，可以发现：在复合圆柱壳制备过程中存在的折返层缺陷对其结构响应及稳定性影响较大，因此在制备过程中需要尽可能提高工艺，保证铜线折返层的完整性及均匀性；通过对0°、10°和15°等3 种螺旋角度复合圆柱壳模型的计算，可以发现复合圆柱壳的螺旋角度对其稳定性影响不大；而由于铜线直径会直接影响复合圆柱壳各层材料的厚度，因此对稳定性影响较大，是进行装置设计及实验需重点考虑的参数。