基于AUTODYN的火工品爆炸冲击响应仿真分析

2018-05-17杨泽雨

导弹与航天运载技术 2018年2期

薛锋，张刚，王飞，韩铭，杨泽雨

0 引言

航天器在飞行过程中，火工品装置爆炸会在邻近结构上产生剧烈的瞬态机械响应[1～4]。而火工品爆炸会在结构上产生剧烈的瞬态冲击加速度响应。这种强烈的高频冲击所产生的高冲击、宽频带的爆炸冲击环境极易造成航天器结构和仪器设备的损伤及故障，从而导致设备无法正常工作，甚至造成发射失败。因此，爆炸冲击环境对航天器系统的仪器设备是不可忽略的力学环境，必须进行相应的地面试验验证，考核仪器设备的工作稳定性。但由于航天器结构尺寸复杂且比较大，每个部段都进行真实试验不仅费时、费力，而且成本极大。目前通常采用火工式爆炸冲击模拟试验来考核产品，从而降低试验成本及周期。

由于火工品爆炸冲击是一种复杂的非线性瞬态响应过程，伴随着高温、高压的物理过程，同时还涉及到炸药与空气、炸药与结构及结构与空气的动态相互耦合作用[4]，因此相关的研究、文献资料少。目前研究主要通过试验提出一些经验规律及通过数值计算软件进行模拟研究，但对于爆炸载荷作用下结构的高频响应计算比较困难。

AUTODYN软件是一个显式非线性动力分析软件，可以解决固体、流体、和液体的动态特性及他们之间的相互作用的高度非线性动力学问题[6]。在航天、航空及军工行业中具有重要的应用价值。爆炸冲击是一个高度非线性的动态问题，属于大变形、强非线性问题，同时还涉及到流体与固体的瞬态耦合作用。在火工品爆炸冲击中，流体材料（空气、炸药）采用欧拉网格，其它结构（响应板）采用拉格朗日网格，然后通过通用接触设置它们之间流固耦合相互作用。本文以AUTODYN软件为计算平台，对火药激励响应板爆炸冲击模拟试验装置的冲击加载过程及响应进行了仿真分析，并得到了响应板上近场、中场、远场的变化规律，进行了响应的试验验证。通过仿真计算分析及相关试验研究，最终实现了高量级火工品爆炸冲击环境。为火工品激励响应板爆炸冲击模拟试验提供了一定的指导。

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（三氟氯乙烯/偏氟乙烯共聚物）。炸药、空气采用3D多物质Euler算法，空气外围施加压力流出边界条件来模拟无限空气域。响应板采用Lagrange求解器，结合火药激励响应板爆炸冲击模拟试验装置实际试验情况定义响应板沿宽度方向固支。模型采用欧拉-拉格朗日耦合作用进行计算，这里选用完全耦合进行计算。

2 材料计算模型参数

2.1 PBX9010炸药的状态方程

爆轰产物中的压力用JWL状态方程[8]来描述：

式中 E为单位质量内能；υ为比容；A，B，R1，R2，ω为常数。其中，式（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

2.2 响应板材料模型

响应板的材料为08F优质碳素结构钢。在本文中，08F优质碳素结构钢的强度模型采用 Johnson-Cook。该强度模型可以考虑高温高压下材料应变硬化、应变率硬化、温度软化的影响，适用于描述金属材料从低应变率到高应变率下的动态行为[8]。

Johnson-Cook强度模型表达式为

式中 C为屈服强度；D，n为应变的强化系数；E1为应变率的敏感系数；npε为等效塑性应变；*pε为塑性应变率；m为温度软化参数。

2.3 空气材料模型

空气采用理想气体模型，模型如下：

式中 p为气体压力；pA为外部环境压力；ρ为初始空气密度；R为气体常数；θ为初始气体温度；θZ为绝对零度。计算中，初始空气域初始温度为300 K，密度为1.177 kg/m3,气体常数为 287 J/（kg•K），外界空气压力pA为 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

图4为响应板在t=1.998 ms时的火药激励位置附近的变形情况与应力云图，在爆炸冲击作用下响应板产生了明显的塑性变形。在该时刻属于爆轰波高压阶段，高频冲击较大，炸药起爆后响应板在 t=1.998 ms时，响应板就已经产生了明显的变形。说明炸药爆炸瞬间，爆炸物质在有限的空间和极短的时间内释放出大量能量，产生的大量高温、高压产物迅速挤压周围介质，使介质的压强、密度、温度等物理量迅速上升从而形成爆炸冲击波，爆炸冲击波与响应板发生相互作用，使响应板发生塑性变形。图5为在t=1.998 ms时的仿真模型变化情况。

图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编程进行了处理，得到了相应的时域图及冲击响应谱曲线图。G1～G4的加速度时程曲线如图 6所示。由图 6可知，随着测点离炸药激励位置距离的增大，响应板加速度响应幅值整体呈递减趋势。在炸药爆炸瞬间，响应板上瞬间激起一高频、高量级冲击信号，随着时间的增加，该信号瞬间衰减。

图6 仿真计算数据处理结果Fig.6 The Data Processing Results of Simulation

续图6

3.2 仿真计算结果对比分析

针对仿真计算结果，采用Matlab编程分别对近区、中区的仿真计算结果进行了冲击响应谱曲线对比分析，结果如图7～ 9所示。

图7 测点G1、G2、G3、G4冲击响应谱曲线对比（近区与中区）Fig.7 The Curve Comparison of the Response Spectrum for the G1、G2、G3、G4 Points(near-field and middle field)

图8 测点G5、G6冲击响应谱曲线对比（近区与中区）Fig.8 The Curve Comparison of the Response Spectrum for the G5、G6 Points (near-field and middle field)

图9 测点G1、G5、G7冲击响应谱曲线对比（近区）Fig.9 The Curve Comparison of the Response Spectrum for the G1、G5、G7 Points（near-field）

由图7～9可知，近区主要由应力波传播支配，其高频响应量级高，冲击响应谱拐点频率在2000～3000 Hz之间，最大值约为25 000 g左右；随着测点离炸药激励位置距离的增大，冲击响应谱整体呈衰减趋势，且其拐点频率趋势变化很小（1000 Hz～3000 Hz之间），低频部分由100 g左右减小至30 g，在进行试验时可通过计算结果选择合适的激励点，避免低频响应对产品造成损坏。随着距离的增加，低频减小较明显，冲击响应谱最大值呈减小趋势。由图9可知，近区的G1、G5、G7位置略有不同，但近区的高低频量级、冲击响应谱趋势，拐点频率一致，从一个方面验证了仿真计算结果的准确性。

可根据仿真计算结果。选择合适的火药激励位置、火药药量等来实现不同的试验要求，减少试验调试时间，减低试验成本。

4 仿真结果与试验结果对比分析

4.1 试验设计

设计了一套火药激励响应板爆炸冲击模拟试验装置（图1）。该装置是以柔性炸药条激励响应板来模拟爆炸冲击环境的试验系统。被试产品安装在响应板的合适位置，在响应板背面安置了柔性炸药条，产品的输入点上利用加速度传感器监测载荷条件。通过调节火药激励位置、形状、波形发生器、响应板边界等来实现不同的冲击响应谱，满足不同的试验要求。

4.2 仿真计算与试验结果对比分析

如上所述，对于高压、高频瞬态爆炸与结构动态耦合作用的计算比较困难，对结构爆炸冲击响应高精度计算就更加困难。因此，根据计算进行了相应的试验研究，应用高频响、高量级冲击传感器对响应板上特征点响应过程进行记录，冲击传感器方向为轴向，并将试验结果与数值计算结果进行了对比，结果如下图10示。图10为某产品的爆炸冲击模拟试验结果，该试验的激励点位置与仿真计算中的G7一致，试验中火药形状为圆环形，导爆索使用80RDXXHV，药量为10.2 g，仿真计算时的火药等效为圆柱状，药量与试验使用相同。仿真计算结果与试验结果的冲击响应谱型整体趋势基本一致，最高量级约为20 000 g，拐点频率在1000 Hz～2000 Hz之间和4000 Hz左右，验证了仿真计算结果的正确性。

通过仿真计算可以得到不同火药形状、药量，激励位置，不同区域的冲击响应谱及拐点频率，在进行试验时可根据不同的试验要求，选择合适的火药激励位置、火药形状及药量、响应板边界等实现爆炸冲击模拟试验。