弹丸刚体外弹道及其引信力学环境数值解算稳定性

引信外弹道环境,不但与弹丸外弹道质心运动有关,还与弹丸外弹道绕心运动有关。计算机数值解算能力和刚体动力学理论的发展,推动了弹丸作为刚体在空中运动的研究,弹道方程也逐渐从质点弹道方程发展为四自由度、五自由度和六自由度刚体外弹道方程,同时也对弹丸的非线性运动赋予较为准确的数学表达,即非线性微分方程组。

无控弹丸外弹道方程组一般是一阶变系数联立方程组,通常只能用数值方法求得数值解,仅在一些特定条件下经过适当的简化才能求得近似解析解[1]。弹道数值解法能以较高的准确度对弹丸外弹道微分方程组进行积分求解,但数值解法相对繁杂[2]。弹道学研究的各个领域以及求解弹道问题的各种方法均有自己特殊的误差,包括近似计算误差、数学误差和取整误差等,其中数学误差与求解方法有关,不同的数值解算方法所产生的误差不同,且不可避免。在数值解算时,积分步长的选择是整个计算中最重要的一步,因为步长的大小不仅决定了计算精度,而且决定了计算工作量[3]。文献[1]提到专业人员常根据计算经验选取步长,如应用质点外弹道方程解算外弹道时,可取步长0.1～0.3 s,对于刚体弹道方程,时间积分步长必须小于0.005 s。文献[2]认为龙格-库塔法的截断误差正比于积分时间步长的5次方,在解算弹道时积分时间步长越小,精度越高,但过小会增加计算时间,并未给出针对龙格-库塔法求解弹道方程的具体积分时间步长。文献[4]认为在考虑弹丸摆动时,求解外弹道方程的积分时间步长应小于1/10的摆动周期。文献[5]认为应用四阶龙格-库塔法对刚体弹道模型和简化刚体弹道模型的外弹道微分方程组解算时,需要较小的积分时间步长,如果积分步长过大,计算结果会发散。文献[6]认为用定步长四阶龙格-库塔法求解刚体外弹道微分方程组的积分时间步长通常为0.005 s,计算简化刚体弹道模型的外弹道微分方程组的积分时间步长通常取为0.5 s。文献[7]认为在同样积分时间步长的条件下,龙格-库塔法的计算精度要比欧拉法高,因此采用龙格-库塔法计算弹丸六自由度外弹道微分方程组获取弹道诸元,但未提及使求解稳定的计算步长。文献[8-10]仅利用刚体弹道方程进行外弹道求解,并未涉及数值解算方法和求解稳定性。目前,尚未见有针对积分时间步长对弹丸刚体弹道方程组解算结果影响的详细仿真研究。

无控弹丸六自由度刚体外弹道常微分方程组简称弹丸刚体外弹道方程,可用作各种工程条件下的弹道计算模型,也可作为射表计算模型,其计算结果最全面,计算精度最高[11]。在利用弹丸刚体外弹道方程解算外弹道绕心运动规律并进一步获取引信外弹道力学环境时,外弹道解算结果和引信惯性零部件轴向过载峰值的稳定性受数值解算方法和积分时间步长影响较大。通过对文献[12]建立的弹丸刚体外弹道模型解算得到的弹道诸元和绕心运动规律以及文献[13]建立的引信惯性零部件轴向过载数学模型求解得到的引信惯性零部件轴向过载峰值进行计算方法研究,对比不同数值解算方法和不同积分时间步长对解算结果的影响,从而找出适用于弹丸六自由度弹道方程和外弹道引信惯性零部件轴向过载数学模型数值解算方法以及经济、准确的积分时间步长。

1 常微分方程组数值解算方法

1.1 龙格-库塔法

龙格-库塔法是以函数y(x)的泰勒级数为基础的改进方法,通过在(xn,xn+1)内预报多个点的斜率值,将其加权平均得到平均斜率的近似值,进而获得更高精度的常微分方程组的数值解,常用为四阶龙格-库塔法。龙格-库塔法的具体表达式见文献[1]。

1.2 阿当姆斯预报-校正法

阿当姆斯预报-校正法属于多步法,其思想为:在已知yn、yn-1、yn-2、yn-3时,充分利用前面多步的信息来预测yn+1,以获得较高的预测精度,其优点在于计算量小,计算速度快,缺点在于不能自启动。阿当姆斯预报-校正法的具体表达式见文献[1]。

1.3 吉尔法

吉尔法是一种向后差分的数值积分方法,适用于求解刚性微分方程组,其优点为:容易改变阶和步长,可应用于高阶和高稳定的格式,解隐式方程组,所需工作量小。文献[14]分别采用吉尔法和龙格-库塔法对引信球转子转正运动常微分方程组进行数值解算,对比发现吉尔法的求解精度和计算速度均优于龙格-库塔法。吉尔法的具体表达式见文献[14]。

2 刚体外弹道方程求解稳定性分析

2.1 不同数值解算方法对弹道解算结果的影响

为研究不同数值解算方法对弹丸六自由度刚体外弹道微分方程组数值解算的影响,应用文献[12]介绍的弹丸六自由度刚体外弹道方程及其应用求解程序(文献[12]已验证该数学模型及其应用求解程序的可信性),分别采用龙格-库塔法、吉尔法和阿当姆斯预报-校正法进行弹丸外弹道解算。仿真所用计算机CPU型号为Intel(R) Xeon(R) Gold 6434 3.7GHz (64位),以某155 mm口径火炮榴弹为研究对象,进行仿真研究。弹丸结构参数均与射表条件一致。发射参数为初速930 m/s、射角49.812°,初始章动角取0°,积分时间步长取0.005 s[1]。应用不同数值解算方法所得到的弹道诸元与射表对比,如表1所示。

由表1可知,无论是四阶龙格-库塔法还是吉尔法或者阿当姆斯预报-校正法,在积分步长为0.005 s时,对于弹丸六自由刚体外弹道方程所解算出的弹道诸元完全一致,与射表之间的误差均在4.4%以下,因射表数据均由多次靶试试验拟合而成,而数值解算用弹丸气动参数均由流体动力学仿真软件Fluent仿真而来[15],两者结果存在偏差无法避免。在解算结果相同的条件下,阿当姆斯预报-校正法解算时间最快,其次是龙格-库塔法,最慢的是吉尔法。但大约均在10 s以内,差别并无实用意义。

2.2 时间步长对弹道诸元的影响

为了获取刚体外弹道方程数值解算对积分时间步长的需求,选用龙格-库塔法针对不同积分时间步长进行仿真研究,仿真参数与表1相同。对比解算后的诸元,即飞行时间、落速、射程、落角、最大高度和偏流的收敛情况,以及达到收敛时所需的最小积分时间步长。弹道诸元收敛情况如图1所示,仿真结果如表2所示。

由图1和表2可知,积分时间步长取0.05 s时,在弹丸六自由度刚体外弹道方程解算得到的弹道诸元中,飞行时间、落速、落角、最大高度和偏流的解算结果均收敛至小数点后第二位,射程的解算结果收敛至个位,达到米级精度;积分时间步长取0.002 5 s时,射程解算结果达到分米级精度;积分时间步长取0.000 5 s时,射程解算结果达到厘米级精度。因此可以认定为:积分时间步长为0.05 s时解算结果已经收敛,此时计算机时在6.77 s左右,对于无控火炮弹丸落点精度要求,若再继续减小积分时间步长对于解算结果意义不大。

通过分析计算机时可知,当积分时间步长较大时(大于0.001 s),计算机时变化不大,决定计算机时的主要为计算机配置和软件程序等客观因素;当积分时间步长较小时(小于等于0.001 s),决定计算机时的主要因素为积分时间步长,并与积分时间步长成负相关关系。

为验证仿真的普适性,增加初速335 m/s、射角30°、其他参数与表1相同的算例,采用龙格-库塔法在不同积分时间步长的情况下计算的弹道诸元,如表3所示。

对比表2和表3可知,不同仿真参数下,弹道诸元收敛情况一致,证明了方法的普适性。在外弹道过程中,质点运动确定后,弹丸整体运动随即确定。因此,在后续研究时间步长对绕心运动和引信力学环境时,仿真结果具有普适性。

2.3 时间步长对绕心运动诸元的影响

弹丸在外弹道飞行过程中的绕心运动诸元可由欧拉角表示,可通过弹丸六自由度刚体外弹道方程解算结果转换得到[3,15],对解算后的章动角和进动角再求1阶导数和2阶导数,就可得到章动角速度、进动角速度、章动角加速度和进动角加速度。通过改变时间步长对外弹道绕心运动进行解算,在整个外弹道过程中选取4个固定的时间点,对比不同积分时间步长时解算结果差异,进而判定解算结果是否收敛。仿真参数与表1的相同。由于弹丸六自由度微分方程组在解算弹丸绕心运动时刚性较强,不得不选用复杂的吉尔方法。弹道上的4个固定观测点分别为弹丸飞出炮口第0.5、5、10、50 s。所得结果如表4～表7所示。其中解算结果数值精度为小数点后保留4位。

由表4—表7可知,章动角、章动角速度和进动角在积分时间步长为0.1 s时已收敛。章动角加速度、进动角速度、进动角加速度在积分时间步长为10-3 s时与10-4 s时的变化量小于5.3%;积分时间步长为10-4 s时与10-5 s时的变化量小于1%,可近似认为在积分时间步长不大于10-4 s时解算结果基本收敛。若要获得较高精度,需选取积分时间步长不大于10-5 s,即10 μs。

3 引信惯性件轴向过载求解稳定性

3.1 时间步长对引信惯性零部件轴向过载异变峰值的影响分析

根据引信惯性零部件轴向过载数学模型[13],轴向过载的理论计算值与外弹道飞行过程中的绕心运动诸元密切相关,因此,以弹丸六自由度刚体外弹道方程解算后的绕心运动诸元为基础计算的引信惯性零部件轴向过载的稳定性也会受积分时间步长的影响。引起引信弹道炸的原因有可能来自于引信惯性零部件轴向过载突变峰值[13],该峰值的解算准确性对引信安全性设计和分析具有重要意义,因此需对不同积分时间步长下的引信惯性零部件轴向过载峰值进行仿真研究。对比3种数值计算方法仿真结果中的首个突变峰值,如表8所示,3种数值解算方法在不同积分时间步长的轴向过载峰值出现时刻和峰值大小收敛情况如图2和图3所示。

由表8、图2和图3可知,龙格-库塔法、吉尔法和阿当姆斯预报-校正法解算后的引信惯性零部件轴向过载峰值及其出现时刻稳定所需的积分时间步长不同。在积分时间步长不大于10-4 s时,龙格-库塔法解算结果中的轴向过载出现时刻已经收敛,吉尔法收敛所需最大积分时间步长为10-5 s,阿当姆斯预报校正法为10-3 s。对于轴向过载峰值的仿真,在积分时间步长不大于10-5 s时,龙格-库塔法和吉尔法均已收敛,阿当姆斯预报-校正法则需不大于10-6 s。在相同条件下,积分时间步长不小于10-5 s时,3种数值解算方法所需计算用时相差不大;积分时间步长小于10-5 s时,积分时间步长越小,吉尔法计算时间越快,其次为龙格-库塔法,阿当姆斯预报-校正法最慢。

从表8看,3种解算方法计算外弹道得出的引信惯性零部件轴向过载首个峰值及其出现时刻均不相同,说明在峰值出现附近,引信惯性零部件轴向过载数学模型存在奇异,此时方程刚性较强。而在这3种解算方法中,只有吉尔法可以求解刚性微分方程,因此龙格-库塔法和阿当姆斯预报-校正法此时及以后会求解失真。为验证此观点,取弹丸飞出炮口第1 s、第2 s和第3 s的过载值进行比较,结果如表9所示。

由表9可知,3种解算方法得到的首个轴向过载峰值最早出现在弹丸飞出炮口后第5.12 s,弹丸飞出炮口1 s、2 s和3 s均在外弹道引信惯性零部件轴向过载突变峰值出现之前。由表9可知,3种解算方法在这3个时刻的过载值相差不大,在数值上的区别为小数点后第4位,彼此相差小于2.2%,说明在突变峰值发生之前,引信惯性零部件轴向过载解算结果较为稳定。

3.2 数值解算方法对轴向过载峰值仿真准确定性分析

为对比不同数值解算方法对引信惯性零部件轴向过载异变峰值仿真准确性,参考文献[13]中靶场射击试验引信活机体部位轴向过载实测值中的首个突变峰值和峰值出现的时刻,对比3种解算方法仿真结果与实测值的误差。不同数值解算方法与实测值对比结果如表10所示,仿真积分时间步长取10-5 s。

由表10可知,3种数值解算方法对引信惯性零部件轴向过载异变峰值的仿真结果中,只有吉尔法的接近实测值,而龙格-库塔法和阿当姆斯预报-校正法仿真结果与实测值差异都较大。

4 结论

1) 在利用弹丸刚体外弹道方程求解外弹道诸元时,龙格-库塔法、吉尔法和阿当姆斯预报-校正法的精度一致,保证求解稳定的最大积分时间步长为0.05 s。

2) 在求解弹丸外弹道绕心运动诸元时,保证解算结果稳定的最大积分时间步长为10-4 s。

3) 在求解外弹道引信惯性零部件轴向过载峰值时,发现龙格-库塔法和阿当姆斯预报-校正法在轴向过载突变发生时刻及以后会求解失真,但吉尔法未失真。

4) 在求解外弹道引信惯性零部件轴向过载峰值时,保证解算结果稳定的最大积分时间步长为10-5 s。

[1] 韩子鹏,常思江,赵子华,等.弹箭外弹道学[M].北京:北京理工大学出版社,2022. HAN Zipeng,CHANG Sijiang,ZHAO Zihua,et al.Exterior ballistics of projectiles and rockets[M].Beijing:Beijing Institute of Technology Press,2022.

[2] 王中元,史金光,常思江,等.现代外弹道学[M].北京:科学出版社,2024. WANG Zhongyuan,SHI Jinguang,CHANG Sijiang,et al.Modern external ballistics[M].Beijing:Science Press,2024.

[3] A.A.德米特里耶夫斯基,Л.H.雷申科,C.C.波哥吉斯托夫.外弹道学[M].韩子鹏.北京:国防工业出版社,2000. DEMITLIYEFU A A,LESTENKO Л H,BOGOZYSTOV C C.Translated by HAN Z P.Exterior ballistics[M].Beijing:National Defense Industry Press,2000.

[4] 浦发.外弹道学[M].北京:国防工业出版社,1980. PU Fa.Exterior ballistics[M].Beijing:National Defense Industry Press,1980.

[5] 王兆胜.火炮射击精度分析的模型与应用[M].北京:国防工业出版社,2013. WANG Zhaosheng.Analysis models of gun fring accuracy and their applications[M].Beijing:National Defense Industry Press,2013.

[6] 郭锡福.远程火炮武器系统射击精度分析[M].北京:国防工业出版社,2004. GUO Xifu.Firing accuracy analysis for long range gun weapon systems[M].Beijing:National Defense Industry Press,2004.

[7] 徐超,王治霖,王武刚,等.基于PSO优化的激光末制导炮弹诸元解算方法[J].弹箭与制导学报,2022,42(2):122-128. XU Chao,WANG Zhilin,WANG Wugang,et al.Calculation method of laser terminal guidance projectile firing data based on PSO[J].Journal of Projectiles,Rockets,Missiles and Guidance,2022,42(2):122-128.

[8] HENDRIK R,MARKUS G,RAINER B.Thermal weapon sights with integrated fire control computers:algorithms and experiences[C]//Infrared Technology and Applications XXXIV pt.1.:SPIE-The International Society for Optical Engineering,2008.

[9] 张百川,毕文豪,张安,等.拓展高原环境的航空炸弹弹道快速解算方法[J].系统工程与电子技术,2024,46(6):2073-2081. ZHANG Baichuan,BI Wenhao,ZHANG An,et al.Fast calculation method of aviation bomb trajectory adapted to plateau environment[J].Systems Engineering and Electronics,2024,46(6):2073-2081.

[10] 李阳,牛志一,刘晓光,等.周向敏感子弹动力学建模及外弹道特性分析[J].兵器装备工程学报,2023,44(6):94-102. LI Yang,NIU Zhiyi,LIU Xiaoguang,et al.Dynamic modeling and exterior trajectory characteristics analysis of transversesensitive submunition[J].Journal of Ordnance Equipment Engineering,2023,44(6):94-102.

[11] 刘世平,管军,常思江,等.弹箭飞行试验与气动参数辨识[M].北京:国防工业出版社,2023. LIU Shiping,GUAN Jun,CHANG Sijiang,et al.Flight test and aerodynamic coefficients identification of projectiles[M].Beijing:National Defense Industry Press,2023.

[12] 彭启蒙,王雨时,项帆,等.基于刚体外弹道学仿真的中大口径旋转炮弹引信力学环境[J].探测与控制学报,2023,45(2):37-44. PENG Qimeng,WANG Yushi,XIANG Fan,et al.Mechanical environment of fuze for medium and large caliber rotating projectile based on rigid body exterior ballistics simulation[J].Journal of Detection &Control,2023,45(2):37-44.

[13] 回建豪.旋转弹丸引信外弹道力学分析与弹道炸影响因素研究[D].南京:南京理工大学,2024. HUI Jianhao.Analysis of external ballistic mechanics of rotating projectile fuze and study on factors influencing ballistic burst[D].Nanjing:Nanjing University of Science and Technology,2024.

[14] 陈会光,王雨时,闻泉.用Gear法求解病态的引信球转子运动常微分方程[J].探测与控制学报,2008,30(4):30-35,39. CHEN Huiguang,WANG Yushi,WEN Quan.Solution to stiff ordinary differential equations on the movement of fuze ball rotor with gear method[J].Journal of Detection &Control,2008,30(4):30-35,39.

[15] 彭启蒙.旋转炮弹引信弹道环境及弹道炸成因研究[D].南京:南京理工大学,2022. PENG Qimeng.Research on The ballistic environment and causes of ballistic explosion of rotating projectiles fuze[D].Nanjing:Nanjing University of Science and Technology,2022.