0 引言
固体火箭发动机作为导弹和运载火箭的重要构件,其内部的高能固体推进剂组分中含有大量高能炸药,在外界作用下极易发生爆炸或爆轰[1],导致殉爆现象的发生。
目前,对于固体火箭发动机和其他含能材料的殉爆问题,许多学者开展了研究。路胜卓等[2]对带壳高能固体推进剂的殉爆过程进行了数值模拟和实验研究,分析了推进剂的殉爆特性。Pang等[3]采用任意拉格朗日欧拉(ALE)方法模拟了纤维复合材料壳体固体火箭发动机的殉爆过程,发现发动机近距离殉爆的主要原因是爆轰产物,复合材料壳体破片的撞击作用有限。陈朗等[4-5]进行了壳装固黑铝炸药的殉爆实验和数值模拟计算,考虑了主发炸药爆炸冲击波对被发炸药的冲击起爆过程。Dong等[6]通过数值模拟方法研究了引信传爆序列的殉爆行为,提出了预测临界殉爆距离的理论模型。
由于开展固体火箭发动机殉爆试验的风险和成本较高,且可重复性差,数值模拟方法是研究其殉爆过程的一种有效手段。现有文献中涉及殉爆的数值模拟大多基于传统的ALE方法,存在模型复杂、网格扭曲和计算效率低的问题。目前,许多学者已将光滑粒子流体动力学(SPH)方法应用于爆炸冲击响应研究中[7-10]。鉴于此,本研究首先分别采用ALE方法和SPH方法对近场殉爆模型进行数值模拟,验证SPH方法计算结果的准确性和计算的高效性。其次,采用SPH方法建立固体火箭发动机的殉爆数值模型,探究主、被发发动机的间隔距离、尺寸效应以及排列方式对其殉爆响应特性的影响。
1 SPH方法
1.1 SPH方法基本理论
SPH方法是一种基于无网格算法的连续介质动力学计算方法,该方法将连续的流体或固体用相互作用的质点组描述,对负载有力学量的流动的粒子进行跟踪,进而求得整个系统的力学行为[11-12]。其本质上是模拟流体流动的一种拉格朗日型粒子方法,因此,其控制方程采用Lagrange形式的Navier-Stockes方程,可通过三大守恒定律得到[13]:
(1)
式(1)中: α, β=1,2,3,代表坐标方向,遵循张量运算字母指标法; ρ、e、vα、vβ、xα、xβ、σαβ、t分别表示密度、内能、速度分量、空间坐标、总应力张量和时间;Fα表示作用在单位质量流体上的体力。
通过SPH的核近似和粒子近似原理,对式(1)进行粒子近似,得到用于爆炸数值模拟的离散化SPH控制方程[14]:
(2)
1.2 SPH方法有效性验证
为证明采用SPH方法模拟殉爆问题的有效性,建立近场殉爆数值模型如图1所示。分别采用ALE方法和SPH方法计算同样的模型,对比分析隔板的动响应和被发装药内部高斯点的压力变化。ALE模型中需建立空气域,主发装药为欧拉网格,隔板和被发装药为拉格朗日网格;SPH模型中,将主、被发装药均离散为粒子,隔板为拉格朗日网格。设置主发装药左侧下端点处起爆,在被发装药中心处设置高斯点,以输出被发装药内部的压力变化曲线。主、被发装药均为COMP B炸药,隔板材料为STEEL 1006钢,材料参数取自AUTODYN标准材料库。
图1 近场殉爆模型
Fig.1 Near-field sympathetic detonation model
图2为分别采用2种方法计算得到的被发装药在不同时刻的压力云图。SPH方法和ALE方法计算得到的被发装药在爆炸冲击作用下的响应过程基本吻合。
图2 ALE方法(左)和SPH方法(右)不同时刻被发装药压力云图
Fig.2 ALE method (left) and SPH method (right) clouds of acceptor charge pressure at different moments
采用SPH方法与ALE方法计算得到的隔板动能-时间历程曲线和被发装药高斯点处压力-时间历程曲线对比如图3所示。2种方法计算得到的隔板动能变化趋势和被发装药高斯点处的压力脉冲持续时间和变化趋势基本吻合。隔板动能第一次上升是由于主发装药爆轰产物的冲击作用,当冲击作用减弱后动能下降,动能再次上升是由于被发装药发生殉爆,被发装药的爆轰产物对隔板产生更加剧烈的冲击作用。
图3 2种方法的隔板动能变化和高斯点压力变化
Fig.3 Change in kinetic energy of the spacer and change in Gaussian point pressure for the two methods
表1为2种方法仿真结果的对比分析,观察表1可知,相较于ALE方法,基于SPH方法建立的殉爆模型网格单元数减少约46.7%,计算时长减少约37.5%,2种方法计算得到的被发装药高斯点最大压力值相对误差约1.9%。
表1 2种方法仿真结果对比
Table 1 Comparison of simulation results between the two methods
仿真方法被发装药响应等级高斯点最大压力值Pm/GPa网格单元数量计算时长/sALE爆轰24.2238900240SPH爆轰23.7520707150
在进行爆炸冲击作用下的结构和装药响应等非接触爆炸仿真分析时,SPH方法相较于传统的ALE方法有较大的优势。由于主发装药和被发装药之间有一段距离,为了模拟主发装药爆炸后爆轰产物在介质中的传播,ALE方法需要设置较大的空气域,而SPH方法的粒子可以自由膨胀,不需要在可能膨胀的区域预先定义粒子,对于相同的物理模型,会节约更多的计算资源,更好地处理爆轰产物的动态变化和非线性效应。同时,在SPH方法中,由于粒子之间不存在网格关系,因此可以避免大变形时网格畸变导致计算步长过小并影响计算精度的问题。
在模拟殉爆问题时,SPH方法不依赖于网格结构,计算效率高,计算精度较好,处理复杂几何形状灵活,适用于本文中爆炸冲击作用下固体火箭发动机的殉爆响应研究。
2 数值模拟
2.1 材料模型
固体火箭发动机壳体和绝热层材料的状态方程、强度模型和侵蚀准则如表2所示。
表2 发动机壳体和绝热层材料模型
Table 2 Model of the motor casing and insulation material
部件材料状态方程强度模型侵蚀准则壳体钢冲击状态方程詹森-库克模型几何应变绝热层橡胶超弹性状态方程超弹性模型几何应变
采用JWL状态方程描述推进剂的未反应物和反应产物[15],状态方程:
(3)
式(3)中:P为爆轰产物的压力,Pa;V为爆轰产物的相对比容;E为爆轰产物的比内能;ω、R1、R2、A和B为推进剂常数。
推进剂冲击起爆特性采用Lee-Tarver三项式点火增长模型描述[16-17],状态方程:
G2(1-λ)eλgpz
(4)
式(4)中: λ为推进剂反应度;t为时间,μs;ρ为密度,kg·m3;ρ0为初始密度,kg·m3;p为压力,GPa;I、b、x为点火控制参量,a为点火临界压缩度参数,G1、c、d、y为点火后早期增长反应控制参数,G2、e、g、z为高压反应控制参数。
推进剂的JWL状态方程参数和点火增长模型参数如表3和表4所示[18]。
表3 推进剂JWL状态方程参数
Table 3 JWL equation of state parameters for propellant
A/GPaB/GPaR1R2ω反应物909.5962.0551.820.2未反应物7000-1.6741010.8
表4 推进剂点火增长模型参数
Table 4 Ignition growth model parameters of propellant
I/μs-1bxaG1c440.22240.011110.222dyG2egz0.6671.662000.3330.6672
2.2 物理模型
运用AUTODYN软件开展固体火箭发动机的殉爆数值模拟研究,由于发动机结构的对称性,建立其1/2模型如图4所示。
图4 固体火箭发动机有限元模型
Fig.4 Solid rocket motor finite element model
通过Hypermesh软件对发动机模型进行网格划分,生成K文件后导入AUTODYN,建立固体火箭发动机的光滑粒子模型。主/被发发动机的结构和尺寸完全一致,发动机的壳体和绝热层采用拉格朗日法描述,推进剂采用SPH方法描述。经多次数值模拟试验,在兼顾计算精度和计算效率的前提下,选取SPH粒子尺寸为0.15 mm。起爆点设置于主发发动机推进剂左侧中心点处,1—16为选取的被发发动机推进剂内部的高斯点,固体火箭发动机殉爆模型如图5所示。
图5 固体火箭发动机殉爆模型
Fig.5 Solid rocket motor sympathetic detonation model
3 结果分析
对于固体火箭发动机殉爆过程中的冲击起爆问题,理想C-J爆轰理论[19]认为,爆轰传播速度趋向于爆炸物的理想爆速,并以该特征速度稳定传播,此时所达到的爆轰状态称为C-J爆轰状态。爆轰状态下,被发发动机的推进剂处于热化学平衡和热力学平衡状态,其压力是一个稳定的数值,称为C-J爆轰压力。因此,在被发发动机受到冲击后,其内部各单元压力保持稳定时,认为达到C-J爆轰压力,可判定被发发动机发生爆轰反应,主要综合被发发动机的压力云图和内部高斯点的压力变化趋势判断响应情况。
3.1 主发发动机爆炸过程
直径为100 mm的主发发动机在爆炸过程中不同时刻的发动机压力云图如图6所示。主发发动机推进剂在起爆点发生爆炸反应并产生爆轰波,爆轰波由起爆点向外持续稳定传播,主发发动机内部压力不断增大,推进剂爆轰产物粒子持续膨胀。
图6 主发发动机起爆压力云图
Fig.6 Donor motor detonation pressure cloud
3.2 间隔距离影响
分别对直径为100 mm的主/被发发动机在平行排列方式下,L1=50、60、70 mm时被发发动机的响应情况进行仿真分析。
图7为L1=50 mm时被发发动机的殉爆压力云图。t=26.50 μs时,在主发发动机爆轰产物的冲击作用下,被发发动机内部压力增长至22.77 GPa,随着压力不断增加,被发发动机内部形成稳定传播的爆轰波,爆轰波由左向右持续传播,在极短时间内引爆整个被发发动机。
图7 L1=50 mm被发发动机压力云图
Fig.7 Pressure cloud of the acceptor motor at L1=50 mm
图8为3种间距下被发发动机内部高斯点处的压力时间历程曲线。由图8(a)可知,L1=50 mm时,被发发动机内部压力峰值持续稳定,发生爆轰反应;由图8(b)可知,被发发动机内部同样形成稳定传播的爆轰波,发生爆轰反应;由图8(c)可知,被发发动机的内部压力较小,最高压力不超过1.2 GPa,低于临界起爆压力,且被发发动机内部压力呈下降趋势,未发生爆轰反应。
图8 平行排列3种间距下被发发动机压力曲线
Fig.8 Pressure curves for three spacings of the acceptor motor when placed in parallell arrangement
3.3 尺寸效应影响
为分析尺寸效应对固体火箭发动机殉爆响应特性的影响,同时将主、被发发动机直径等比例增加至150 mm,分析平行排列方式下,L1=70、80、100 mm时被发发动机的响应程度。图9显示了L1=70 mm时被发发动机在不同时刻的压力云图。
图9 直径150 mm被发发动机在L2=70 mm时压力云图
Fig.9 Pressure cloud for the 150 mm acceptor motor at L2=70 mm
3种间距下被发发动机内部的压力时间历程曲线如图10所示。直径150 mm的被发发动机在间距L1=70、80 mm下发生爆轰反应,发动机内部形成稳定传播的爆轰波,在间距L1=100 mm时,被发发动机内部最高压力值约为1.2 GPa,且随着冲击波的传递,压力大小呈下降趋势,被发发动机未发生爆轰反应。对比可知,同时增加主、被发发动机的直径,将增大主发发动机的爆轰产物强度,使被发发动机受到冲击的面积和能量增大,进而增加发动机的临界殉爆距离。
图10 直径150 mm被发发动机压力曲线
Fig.10 Pressure curves of the 150 mm acceptor motor
3.4 排列方式影响
为分析主/被发发动机排列方式对其殉爆响应特性的影响,将直径100 mm的主/被发发动机相对垂直排列,L2=50、60、70 mm。
图11为L2=50 mm时被发发动机在不同时刻的殉爆响应压力云图。t=27 μs时,被发发动机左侧端部在主发发动机爆轰产物的冲击作用下发生变形,随着爆轰产物的持续作用,被发发动机内部形成稳定传播的爆轰波,被发发动机发生爆轰反应。图12为3种间距下被发发动机内部高斯点处的压力-时间历程曲线。由图12(a)可知,L2=50 mm时,被发发动机内部形成稳定传播的爆轰波,发生爆轰反应;由图12(b)和(c)可知,L2=60、70 mm时,被发发动机均未发生爆轰反应,内部最高压力远低于爆轰压力。
图11 L2=50 mm被发发动机压力云图
Fig.11 Pressure cloud of the acceptor motor at L2=50 mm
图12 垂直排列3种间距下被发发动机压力曲线
Fig.12 Pressure curves for three spacings of the acceptor motor when placed in vertical arrangement
固体火箭发动机的排列方式对其殉爆响应特性影响显著。由于固体火箭发动机的排列方式不同,主发发动机爆炸产生的爆轰产物对被发发动机的作用效果也不同,平行排列时,主/被发发动机之间的能量传递更为直接,被发发动机受到爆轰产物冲击的面积更大,冲击作用更强,从而增加了临界殉爆距离。而垂直排列时爆轰产物的能量传递路径更加复杂,主发发动机的爆轰产物对被发发动机作用面积较小,降低了临界殉爆距离。
表5为上述不同计算条件下,被发发动机的殉爆响应数值模拟结果。
表5 不同条件下被发发动机响应结果
Table 5 Response results of the acceptor motor under different operating conditions
发动机直径/mm排列方式间隔距离L/mm被发发动机响应程度100平行垂直506070506070爆轰爆轰未爆轰爆轰未爆轰未爆轰150平行7080100爆轰爆轰未爆轰
4 结论
本文采用SPH方法建立了固体火箭发动机的殉爆数值模型,分析了不同条件下被发发动机在主发发动机爆轰产物冲击作用下的响应程度,结论如下:
1) 在开展近场殉爆数值模拟研究时,相较于ALE方法,采用SPH方法在保证较高计算精度的同时,计算效率提高了约37.5%,且不会出现拉格朗日网格畸变等不良情况。
2) 主、被发发动机间距越大,被发发动机受主发发动机爆轰产物的影响程度越低,被发发动机的安全性越高。100 mm直径发动机平行排列,间距50 mm和60 mm时发生殉爆,间距70 mm时未殉爆;100 mm直径发动机垂直排列,间距50 mm时发生殉爆,间距60 mm和70 mm时未殉爆,且间距70 mm时被发发动机的响应程度低于间距60 mm时的响应程度;150 mm直径发动机平行排列,间距70 mm和80 mm时发生殉爆,间距100 mm时未殉爆。
3) 主、被发发动机的尺寸大小对其殉爆响应特性影响较大,同时增加主、被发发动机的直径,将显著提高被发发动机的殉爆响应程度,增加临界殉爆距离。平行排列时,100 mm直径发动机的临界殉爆距离在60～70 mm之间,而150 mm直径发动机的临界殉爆距离在80-100 mm之间。
4) 排列方式对固体火箭发动机的殉爆响应特性有显著影响,相较于垂直排列方式,在平行排列方式下,被发发动机受到主发发动机爆轰产物的冲击作用更强,殉爆响应更剧烈。100 mm直径发动机在平行排列方式下的临界殉爆距离在60～70 mm之间,大于垂直排列方式下的50～60 mm。
[1] 秦能,廖林泉,范红杰,等.几种典型固体推进剂的危险性能实验研究[J].含能材料,2010,18(3):324-329.QIN Neng,LIAO Linquan,FAN Hongjie,et al.Sensitivity performances of several typical solid propellants[J].Chinese Journal of Energetic Materials,2010,18(3):324-329.
[2] 路胜卓,罗卫华,陈卫东,等.壳装高能固体推进剂的殉爆实验与数值模拟[J].哈尔滨工程大学学报,2014,35(12):1507-1511,1552.LU Shengzhuo,LUO Weihua,CHEN Weidong,et al.Experiments and numerical simulations of sympathetic detonation of high-energy solid propellant in shell[J].Journal of Harbin Engineering University,2014,35(12):1507-1511,1552.
[3] PANG S L,CHEN X,XU J S.Numerical simulations of sympathetic detonation of solid rocket motors[J].Journal of Physics:Conference Series,2022,2235(1):012014.
[4] 陈朗,王晨,鲁建英,等.炸药殉爆实验和数值模拟[J].北京理工大学学报,2009,29(6):497-500,524.CHEN Lang,WANG Chen,LU Jianying,et al.Experiment &simulation of sympathetic detonation tests[J].Transactions of Beijing Institute of Technology,2009,29(6):497-500,524.
[5] CHEN L,WANG C,FENG C G,et al.Study on random initiation phenomenon for sympathetic detonation of explosive[J].Defence Technology,2013,9(4):224-228.
[6] DONG L Y,WANG Z J,XIAO Y C,et al.Sympathetic detonation reaction behavior of a fuze explosive train[J].Combustion,Explosion,and Shock Waves,2021,57(5):607-617.
[7] 刘粟涛,周云波,张明,等.爆炸冲击波与破片作用下车辆底部结构动响应数值仿真[J].兵器装备工程学报,2022,43(5):56-62.LIU Sutao,ZHOU Yunbo,ZHANG Ming,et al.Numerical simulation of dynamic response of vehicle bottom structure under action of explosion shock wave and fragments[J].Journal of Ordnance Equipment Engineering,2022,43(5):56-62.
[8] 强洪夫,孙新亚,王广,等.钢箱内部爆炸破坏的SPH数值模拟[J].爆炸与冲击,2019,39(5):24-32.QIANG Hongfu,SUN Xinya,WANG Guang,et al.Numerical simulation on steel box damage under internal explosion by smoothed particle hydrodynamics[J].Explosion and Shock Waves,2019,39(5):24-32.
[9] 耿航.改进SPH方法在冲击载荷下板壳结构动响应中的应用研究[D].哈尔滨:哈尔滨工程大学,2020.GENG Hang.Application of the improved SPH method in the dynamic response of shell structures under impact loads[D].Harbin:Harbin Engineering University,2020.
[10] 梁博,蒋宏业,徐涛龙,等.基于SPH-FEM耦合算法的埋地输气管道近场爆炸冲击动力响应[J].石油学报,2017,38(11):1326-1334.LIANG Bo,JIANG Hongye,XU Taolong,et al.Impact dynamic response of near-field explosion in buried gas pipeline based on SPH-FEM coupling algorithm[J].Acta Petrolei Sinica,2017,38(11):1326-1334.
[11] 杨刚,韩旭,龙述尧.应用SPH方法模拟近水面爆炸[J].工程力学,2008(4):204-208,213.YANG Gang,HAN Xu,LONG Shuyao.Simulation of underwater explosion near air-water surface by SPH method[J].Engineering Mechanics,2008(4):204-208,213.
[12] 张志春,强洪夫,高巍然.一种新型SPH-FEM耦合算法及其在冲击动力学问题中的应用[J].爆炸与冲击,2011,31(3):243-249.ZHANG Zhichun,QIANG Hongfu,GAO Weiran.A new coupled SPH-FEM algorithm and its application to impact dynamics[J].Explosion and Shock Waves,2011,31(3):243-249.
[13] 初文华,明付仁.三维SPH算法在冲击动力学中的应用[M].北京:科学出版社,2017:13-17.CHU Wenhua,MING Furen.Application of the 3D SPH algorithm to impact dynamics[M].Beijing:Science Press,2017:13-17.
[14] XU J X,LIU X L.Analysis of structural response under blast loads using the coupled SPH-FEM approach[J].Journal of Zhejiang University SCIENCE A,2008,9(9):1184-1192.
[15] LEE E,FINGER M,COLLINS W.JWL equation of state coefficients for high explosives[R].Lawrence Livermore National Laboratory,W-7405-ENG-48:1973.
[16] TARVER C M,HALLQUIST J O,ERICKSON L M.Modeling short pulse duration shock initiation of solid explosives[C]//Proceedings of the 8th International Symposium on Detonation.Albuquerque,NM,US:Office of Naval Research,1985:951-961.
[17] 张泽远,沈欣,莫展,等.冲击载荷下固体火箭发动机起爆过程数值模拟[J].弹箭与制导学报,2022,42(3):55-60.ZHANG Zeyuan,SHEN Xin,MO Zhan,et al.Numerical research of shock initiation process of SRM with mechanical impact loading[J].Journal of Projectiles,Rockets,Missiles and Guidance,2022,42(3):55-60.
[18] 崔浩,郭锐,宋浦.固体火箭发动机跌落安全性数值分析[J].兵工学报,2018,39(S1):66-71.CUI Hao,GUO Rui,SONG Pu.Numerical analysis of safety of solid rocket engine during falling process[J].Acta Armamentarii,2018,39(S1):66-71.
[19] 曾祥国.冲击与爆炸动力学基础及应用[M].成都:四川大学出版社,2017:25-29.ZENG Xiangguo.Fundamentals and applications of shock and explosive dynamics[M].Chengdu:Sichuan University Press,2017:25-29.