东方鲎2018年第2期 导 弹 与 航 天 运 载 技 术 No.2 2018 总第360期 MISSILES AND SPACE VEHICLES Sum No.360
收稿日期：2017-01-23；修回日期：2017-09-23
成长的烦恼初中作文作者简介：薛 锋（1985-），女，工程师，主要研究方向为结构动力学、动力学环境试验 文章编号：1004-7182(2018)02-0115-06 DOI ：10.7654/j.issn.1004-7182.20180223
薛 锋，张 刚，王 飞，韩 铭，杨泽雨
（北京强度环境研究所，北京，100076）
摘要：根据火工品爆炸冲击环境模拟装置的基本原理，采用AUTODYN 数值计算软件，对火工品激励响应板爆炸冲击模拟装置进行了仿真研究。得到了不同位置的加速度时程曲线，将计算结果和试验数据进行了对比。结果表明，试验与计算结果吻合较好。通过前期仿真计算减少了试验调试时间与难度，降低了试验成本，提高了火工品爆炸冲击试验的精度，对高量级火工品爆炸冲击试验具有重要指导意义。 关键词：爆炸冲击；火工装置；爆炸冲击试验；AUTODYN 数值模拟；流-固耦合 中图分类号：V416.5 文献标识码：A
Research on the Techniques of Pyrotechnic Explosive Impact Test by
AUTODYN
Xue Feng, Zhang Gang, Wang Fei, Han Ming, Yang Ze-yu
(Beijing Institute of Structure and Environment Engineering, Beijing, 100076)
Abstract: According to the basic theory of pyrotechnic explosive impact test, the simulation for a pyrotechnically excitation
device simulator is carried out by use of the commercial software AUTODYN. Several of the acceleration time history curves of different locations have been simulated, compared the calculated results with experimental data. The results show that simulation and the experimental results are consistent with each other rather well. The experiment cost and debugging time is reduced according to the prophase numerical simulation, and improved the precision control technique, provides effective guide for the pyrotechnic explosive impact test.
Key words: Explosive impact; Pyrotechnics; Explosive impact test; AUTODYN numerical simulation; Fluid-solid coupling
0 引 言
航天器在飞行过程中，火工品装置爆炸会在邻近结构上产生剧烈的瞬态机械响应[1~4]。而火工品爆炸会在结构上产生剧烈的瞬态冲击加速度响应。这种强烈的高频冲击所产生的高冲击、宽频带的爆炸冲击环境极易造成航天器结构和仪器设备的损伤及故障，从而导致设备无法正常工作，甚至造成发射失败。因此，爆炸冲击环境对航天器系统的仪器设备是不可忽略的力学环境，必须进行相应的地面试验验证，考核仪器设备的工作稳定性。但由于航天器结构尺寸复杂且比较大，每个部段都进行真实试验不仅费时、费力，而且成本极大。目前通常采用火工式爆炸冲击模拟试验来考核产品，从而降低试验成本及周期。
由于火工品爆炸冲击是一种复杂的非线性瞬态响应过程，伴随着高温、高压的物理过程，同时还涉及
到与空气、与结构及结构与空气的动态相互耦合作用[4]，因此相关的研究、文献资料少。目前研究主要通过试验提出一些经验规律及通过数值计算软件进行模拟研究，但对于爆炸载荷作用下结构的高频响应计算比较困难。
AUTODYN 软件是一个显式非线性动力分析软件，可以解决固体、流体、和液体的动态特性及他们之间的相互作用的高度非线性动力学问题[6]。在航天、航空及军工行业中具有重要的应用价值。爆炸冲击是一个高度非线性的动态问题，属于大变形、强非线性问题，同时还涉及到流体与固体的瞬态耦合作用。在火工品爆炸冲击中，流体材料（空气、）采用欧拉网格，其它结构（响应板）采用拉格朗日网格，然后通过通用接触设置它们之间流固耦合相互作用。本文以AUTODYN 软件为计算平台，对火药激励响应板
导 弹 与 航 天 运 载 技 术 2018年
116
爆炸冲击模拟试验装置的冲击加载过程及响应进行了仿真分析，并得到了响应板上近场、中场、远场的变化规律，进行了响应的试验验证。通过仿真计算分析及相关试验研究，最终实现了高量级火工品爆炸冲击环境。为火工品激励响应板爆炸冲击模拟试验提供了一定的指导。
1 仿真计算模型
1.1 火药激励响应板模型
仿真计算模型以火药激励响应板爆炸冲击模拟试验装置试验为基础，其中响应板为1.2 m×0.6 m×0.01 m 的钢板，材料选用08F 优质碳素结构钢。根据试验研究和分析，火药选用 PBX9010，质量为10.2 g ，结合试验经计算采用直径为18.18 mm ，厚度为5.5 mm 的圆柱体火药模型。试验装置如图1所示。
图1 火药激励响应板爆炸冲击模拟试验装置
Fig.1 The Pyrotechnically Excitation Device Simulator
1.2 网格密度选取
AUTODYN 软件提供了拉格朗日（Lagrange ）、欧拉（Euler ）、任意拉格朗日欧拉（Arbitrary Lagrangian -Eulerian ，ALE ）、光滑粒子流体动力（Smoothed Part icle Hydrodynamics ，SPH ）
和多物质流固耦合方法等多个求解器[7]。爆炸冲击是一个高度非线性的动态问题，属于大变形、强非线性问题，同时还涉及到空气爆炸载荷、响应板与工装结构等的相互耦合作用。由于多物质流固耦合方法中，一些子域用欧拉法，另一些子域上用拉格朗日法，子域分离或者重叠[8]。在火工品爆炸冲击中，流体材料（空气、）采用欧拉和ALE 算法，对其他的结构（响应板）采用拉格朗日算法，然后通过流固耦合方式来处理相互作用。将流固分开建模，并通过流固耦合的方式来处理相互作用，能方便的建立爆炸模型。这里我们选用多物质流固耦合方法。
火药激励响应板模型如图2所示，响应板采用Lagrange 方法，实际划分的网格尺寸大致为 5 mm×5 mm×2 mm ，单元总数为36 000个，同时在模型中建立了7个高斯观测点进行数据的记录，响应板外围与接触处空气采用Euler 算法，Euler 网格尺寸应小于Lagrange 网格，为了获得准确的耦合相互作用，开展了多种工况的仿真计算分析：Euler 和Lagrange 网格尺寸比为1：1、1：2、1：4、1：5及1：7五种工况，经计算分析发现，当Euler 和Lagrange 网格尺寸比例为1：1至1：4和1：5至1：7之间时，误差较大，而在1：4和1：5之间时误差较小，故Lagrange 网格尺寸应选为4~5倍的Euler 网格尺寸。
计算过程中发现Euler 网格过大，甚至Euler 网格尺寸大于Lagrange 网格尺寸时，
Euler 网格会穿过Lagrange 网格发生耦合泄露，经过网格细化之后这种现象得到了很好的抑制。由于
空气域尺寸大细网格会带来模型网格数量大幅度增加和计算时间步长变短，对于较大的模型，计算所需的时间耗费巨大，所以只对响应板与接触处的空气网格加密，加密计算模型如图3所示。
图2 火药激励响应板模型
Fig.2 The Pyrotechnically Excitation Plate Models
图3 火药激励响应板流固耦合模型（1/4模型）
Fig.3 The Fluid-solid Coupling Models of Pyrotechnically
Excitation Plate
1.3 三维流固耦合模型
本文以火药激励响应板爆炸冲击模拟试验装置作为计算模型，仿真计算用AUTODYN 显式非线性动力分析软件。根据模型的对称性，计算模型只取1/4模型，在该模型中，响应板结构采用体单元而不是壳单元。
仿真计算用为PBX9010，质量为10.2 g ，
PBX9010的配方为90%（质量分数，下同）RDX （）和10%Kel-F3700（三氟氯乙烯/偏氟乙烯共聚
薛 锋等 基于AUTODYN 的火工品爆炸冲击响应仿真分析
117
第2期
物）。、空气采用3D 多物质Euler 算法，空气外围施加压力流出边界条件来模拟无限空气域。响应板采用Lagrange 求解器，结合火药激励响应板爆炸冲击模拟试验装置实际试验情况定义响应板沿宽度方向固支。模型采用欧拉-拉格朗日耦合作用进行计算，这里选用完全耦合进行计算。
2 材料计算模型参数
2.1 PBX9010的状态方程
爆轰产物中的压力用JWL 状态方程[8]来描述：
121211R R E
p A e B e
R R υυωωωυυυ−−⎛⎞⎛⎞=−+−+⎜⎟⎜⎟⎝⎠⎝
⎠ （1） 式中 E 为单位质量内能；υ为比容；A ，B ，R 1，R 2，ω为常数。其中，式（1）右端第1项在高压段起主要作用，
第2项在中压段起主要作用，第3项代表低压段[6,9]。在爆轰产物膨胀的后期，方程式前2项的作用可以忽略。PBX9010爆轰产物JWL 状态方程参数，如表1所示。
表1 PBX9010爆轰产物JWL 状态方程参数
Tab.1 The Parameters of JWL Equation of State Detonation
Product for PBX9010
每天改变一点点参数
A /GPa
B /GPa
1R 2R
ω
数值 581.45 6.801 4.1 1.0 0.35
2.2 响应板材料模型
响应板的材料为08F 优质碳素结构钢。在本文中，08F 优质碳素结构钢的强度模型采用Johnson-Cook 。该强度模型可以考虑高温高压下材料应变硬化、应变率硬化、温度软化的影响，适用于描述金属材料从低应变率到高应变率下的动态行为[8]。
客房管理系统设计
Johnson-Cook 强度模型表达式为
()()()*11ln 1n
m y p p H C D E T σεε=++− （2）
式中 C 为屈服强度；D ，n 为应变的强化系数；E 1为
应变率的敏感系数；n
p ε为等效塑性应变；*p ε为塑性应变率；m 为温度软化参数。
2.3 空气材料模型
空气采用理想气体模型，模型如下：
()Z A p p R ρθθ+=− （3） 式中 p 为气体压力；p A 为外部环境压力；ρ为初始空气密度；R 为气体常数；θ为初始气体温度；θZ 为绝对零度。计算中，初始空气域初始温度为300 K ，密度为1.177 kg/m 3，气体常数为287 J/（kg•K ），外界空气压力
p A 为101 340 Pa 。
3 仿真计算结果分析
根据前期大量经验积累，使用直径为36.36 mm ，厚度为5.5 mm 的圆柱体火药激励响应板中心位置，其中为圆柱形，高5.5 mm ，半径为18.18 mm ，爆心（Detonation ）坐标为（0，0，-5.5）位于响应板中心处，对响应板的瞬态冲击响应过程进行了数值模拟。计算时，在响应板的不同位置设置了监测点，研究响应板的近区、中区、远区等不同位置的爆炸冲击响应及其传递规律，获取冲击环境变化规律，从而指导火工品爆炸冲击试验，优化装置的设计，缩短调试时间， 降低试验成本，实现高量级、高频率火工品爆炸冲击环境的可控性及工程实用性。响应板7个监测点的具体位置见表2。
表2 测点位置区域划分及特点
Tab.2 Region Partition and Characteristics of Measuring Points
测点
坐标
离爆炸源距离/mm
区域
环境特点
冲击加速度幅值/g
频率 范围 kHz
G1(100,15,10)81.82 近区
由应力波传播支配 >5000 >100
G2(200,15,10)181.82中区
由应力波
传播及结构响应 支配
1000~5000 >10
G3(300,15,10)281.82中区 — — — G4(400,15,10)381.82中区 — — — G5(18,100,10)81.82 近区 — — — G6
(18,200,10)
181.82中区 — — — G7(100,100,10)
<150
近区
—
—
—
图4为响应板在t =1.998 ms 时的火药激励位置附
近的变形情况与应力云图，在爆炸冲击作用下响应板产生了明显的塑性变形。在该时刻属于爆轰波高压阶段，高频冲击较大，起爆后响应板在t =1.998 ms 时，响应板就已经产生了明显的变形。说明爆炸瞬间，爆炸物质在有限的空间和极短的时间内释放出大量能量，产生的大量高温、高压产物迅速挤压周围介质，使介质的压强、密度、温度等物理量迅速上升从而形成爆炸冲击波，爆炸冲击波与响应板发生相互作用，使响应板发生塑性变形。图5为在t =1.998 ms 时的仿真模型变化情况。
导
弹 与 航 天
运 载 技 术 2018年
本报内部消息
118
图4 响应板应力云图
Fig.4 The Stress Nephogram of Plate
图5 t =3.563ms 时的仿真模型变化情况
Fig.5 The Changes of Simution Models When t =3.563ms
3.1 仿真计算结果数据处理分析
该模型共设置了7个Gauge 点，每个监测点距离火药激励位置不同，通过计算得到了近场、中场的数据。针对仿真计算数据，通过数据处理软件Matlab 编程进行了处理，得到了相应的时域图及冲击响应谱曲线图。G 1~G 4的加速度时程曲线如图6所示。由图6可知，随着测点离激励位置距离的增大，
响应板加速度响应幅值整体呈递减趋势。在爆炸瞬间，响应板上瞬间激起一高频、高量级冲击信号，随着时间的增加，该信号瞬间衰减。
0.002
0.0040.0060.008
-6000
-4000-2000020004000
6000
80001000012000时间/s
冲击加
速度
/g
a
）测点G 1
图6
仿真计算数据处理结果
Fig.6 The
Data
Processing Results of Simulation
冲
击
加速度/
g
b ）测点
G
2
c ）测点
G
3
d ）测点G 4 续图6
3.2 仿真计算结果对比分析
针对仿真计算结果，采用Matlab 编程分别对近区、中区的仿真计算结果进行了冲击响应谱曲线对比分析，结果如图7~ 9所示。
102
103
1010
2
103104频率/Hz
冲击谱幅值/g
105测点G 1
测点G 2
测点G 3
测点G 4
图7 测点G 1、G 2、G 3、G 4冲击响应谱曲线对比
（近区与中区）
Fig.7 The Curve Comparison of the Response Spectrum for the
广州市交运医院
G 1、G 2、G 3、G 4 Points(near-field and middle field)
薛 锋等 基于AUTODYN 的火工品爆炸冲击响应仿真分析
119
第2期
图8 测点G 5、G 6冲击响应谱曲线对比（近区与中区）
Fig.8 The Curve Comparison of the Response Spectrum for the
G 5、G 6 Points (near-field and middle field)
102
100
1010
2
10
3104响应/Hz
测点G 1
测点G 5
测点G 7
图9 测点G 1、G 5、G 7冲击响应谱曲线对比（近区）
Fig.9 The Curve Comparison of the Response Spectrum for the
G 1、G 5、G 7 Points （near-field ）
由图7~9可知，近区主要由应力波传播支配，其高频响应量级高，冲击响应谱拐点频率在2000~3000 Hz 之间，最大值约为25 000 g 左右；随着测点离激励位置距离的增大，冲击响应谱整体呈衰减趋势，且其拐点频率趋势变化很小（1000 Hz~3000 Hz 之间），低频部分由100 g 左右减小至30 g ，
在进行试验时可通过计算结果选择合适的激励点，避免低频响应对产品造成损坏。随着距离的增加，低频减小较明显，冲击响应谱最大值呈减小趋势。由图9可知，近区的G 1、G 5、G 7位置略有不同，但近区的高低频量级、冲击响应谱趋势，拐点频率一致，从一个方面验证了仿真计算结果的准确性。
可根据仿真计算结果。选择合适的火药激励位置、火药药量等来实现不同的试验要求，减少试验调试时间，减低试验成本。
4 仿真结果与试验结果对比分析
4.1 试验设计
设计了一套火药激励响应板爆炸冲击模拟试验装置（图1）。该装置是以柔性条激励响应板来模拟爆炸冲击环境的试验系统。被试产品安装在响应板的
合适位置，在响应板背面安置了柔性条，产品的输入点上利用加速度传感器监测载荷条件。通过调节火药激励位置、形状、波形发生器、响应板边界等来实现不同的冲击响应谱，满足不同的试验要求。 4.2 仿真计算与试验结果对比分析
如上所述，对于高压、高频瞬态爆炸与结构动态
耦合作用的计算比较困难，对结构爆炸冲击响应高精度计算就更加困难。因此，根据计算进行了相应的试验研究，应用高频响、高量级冲击传感器对响应板上特征点响应过程进行记录，冲击传感器方向为轴向，并将试验结果与数值计算结果进行了对比，结果如下图10示。图10为某产品的爆炸冲击模拟试验结果，该试验的激励点位置与仿真计算中的G 7一致，试验中火药形状为圆环形，使用80RDXXHV ，药量为10.2 g ，仿真计算时的火药等效为圆柱状，药量与试验使用相同。仿真计算结果与试验结果的冲击响应谱型整体趋势基本一致，最高量级约为20 000 g ，拐点频率在
1000 Hz~2000 Hz 之间和4000 Hz 左右，验证了仿
真计算结果的正确性。
通过仿真计算可以得到不同火药形状、药量，激励位置，不同区域的冲击响应谱及拐点频率，在进行试验时可根据不同的试验要求，选择合适的火药激励
位置、火药形状及药量、响应板边界等实现爆炸冲击模拟试验。
10
101010响应/Hz
图10 仿真计算与试验结果响应谱对比
Fig.10 The Response Spectrum Comparison of Simulation and Test
5 结 论
本文利用AUTODYN 接触响应板柔性爆炸冲击响应问题进行了数值模拟及试验研究，得出了以下结论：
a ）计算过程中发现欧拉网格过大，甚至欧拉网格尺寸大于拉格朗日网格尺寸时，欧拉网格会穿过朗格朗日网格发生耦合泄露，Euler 网格尺寸与Lagrange 网格尺寸比例为1：4至1：5之间时，误差最小，过大或者过小都会造成误差过大，经过网格细化之后这种现象得到了很好的抑制。
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议