1 引言
由于火炮在现代战争中占据的重要地位，许多国家都耗费大量的人力物力去研究火炮的内弹道规律[1-5]。中北大学研制出放入式电子测压器来获取膛内的压力曲线，其用膛内某一点的压力变化以近似代替整个膛内的压力变化[6]。但是膛内压力分布并不均匀，某些位置处压力变化差异明显，因此很有必要在膛内寻找多个测压点，以便能够获取完整的膛内压力变化情况[7]。
王亚军、舒跃飞等[8-9]提出将火炮身管近似为规则的柱形密闭容器，在建模时利用二维模型进行仿真，并且选取特殊点观察膛压变化。刘小龙等[10]通过ANSYS软件对火药燃烧进行仿真，观察不同时刻火药气体的压力云图，为火炮膛内的压力分布规律提供了技术支持。他们在建模中都将身管模型简化为二维模型，且所选取的位置并不能完全反应整个膛压的变化规律。本文在不考虑弹丸挤进膛线的情况下，将药室容积简化为密闭容器。利用LS-DYNA软件对火药燃烧进行仿真，并在其内选择多个点观察其压力变化情况。
2 仿真模型
2.1 任意拉格朗日-欧拉算法
在处理大位移、大变形以及大应变等非线性问题时，LS-DYNA更能显示出其软件的优势。其算法包含拉格朗日(lagrange)算法，欧拉(euler)算法以及任意拉格朗日-欧拉ALE(arbitrary lagrangian-eulerian，ALE)算法，其中以拉格朗日算法为主[11]。ALE算法结合了上述2种算法的优点，在间隔一定的时间步后，根据材料物质流动情况,再重新构造一组不同于之前的物质以及空间坐标系的网格，网格点也可以随意移动。这样一来既解决了网格畸变的问题，又可以捕捉到物质的边界，便于边界条件的施加[12-13]。
2.2 仿真模型的建立
本文采用了ALE算法对空气以及炸药进行网格划分。图1为建立火炮身管模型与并未建立火炮身管模型几个相同位置处的压力曲线。图1中,A1—D1为未建立身管模型的压力-时程曲线，A2—D2为建立身管模型的压力-时程曲线。
由图1可以看出，建立身管与否对压力曲线的影响较小，因此在建模时为了简化计算模型，将空气模型边界定义为固定边界，故省去了身管模型的建立。
图2为利用测压器在某火炮内的实测曲线与其仿真结果的对比图，仿真曲线是在未考虑弹丸挤进压力的条件下进行的仿真测试，从吻合度来讲，不考虑挤进压力的情况下，后续研究对实际应用也具有一定的参考价值。
图1 不同模型相同位置处的压力-时程曲线
Fig.1 Pressure-time history curves at the same location for different models
图2 实测结果与仿真结果压力-时程曲线
Fig.2 Comparison of pressure-time history curve between measured results and simulation results
药室参数如表1所示[14]，药筒各部分参数如图3所示，建立了85 mm口径的某火炮药室，并且去掉了药室内的一些特征参数，建立的药室模型如图4所示。
空气用*MAT_NULL材料模型以及线性多样式方程，参数如表2所示。
表1 85 mm口径某火炮药筒结构尺寸
Table 1 Structure and dimension of a certain artillery cartridge with 85 mm caliber
不同口径药筒/mm药筒口部/mm内径d1厚度t1药筒斜肩下部/mm外径d2厚度t2药筒根部/mm外径d3厚度t3药筒底缘/mm外径d4厚度t4药筒底部厚度t5/mm药筒长度l/mm内跟脚半径R/mm8583.201.3094.701.50101.702.85112.004.008.00630.0010.00
图3 药筒各部位尺寸
Fig.3 Dimensions of each part of the cartridge
图4 口径为85 mm某火炮药室模型
Fig.4 Model of an artillery charge room with a caliber of 85 mm
表2 空气参数
Table 2 Air parameters
ρ/(g·cm-3)C4C5E0/JV01.29×10-30.400.402.50×10-61.00
状态方程如下：
P=C0+C1μ+C2μ2+C3μ3+
(C4+C5μ+C6μ2)E
(1)
式(1)中: μ=(ρ/ρ0)-1， ρ为当前密度，ρ0为初始密度;E为炸药的内能；C0、C1、C2、C3、C4、C5、C6为状态方程参数。
炸药采用*MAT_HIGH_EXPLOSIVE_BURN材料模型以及JWL状态方程，参数如表3所示。
状态方程为：
(2)
式(2)中：A、B、R1、R2、ω为JWL状态方程的参数，其值由实验确定；E为炸药的内能；V为当前相对体积。采用等效硝铵炸药装药进行数值仿真。
表3 炸药参数
Table 3 Parameters of explosives
ρ/(g·cm-3)D/(cm·μs-1)PCJ/MPaA/GPaB/GPa1.2000.0690.2703.7103.930R1R2ΩE0/JV04.8200.9700.3500.7501.000
由于并没有在建模的时候就建立身管的模型，因此需要将空气的边界进行固定，炸药是采用物质填充的方法，利用关键字*INITIAL_VOLUME_FRACTION_GEOMETRY来实现。
2.3 仿真结果分析
模型设置的仿真时间为3 000 μs，时间步长为20 μs。图5分别为20、40、60、80、100、120、140、180 μs时刻火药燃烧后的气体压力云图。
图5 不同时刻的火药气体压力云图
Fig.5 Cloud map of gun powder gas pressure at different times
从图5可以看出,在很短的时间内，火药燃烧后的产物要比周围空气的压力高出很多倍，会急剧压缩周围空气，从而使其密度、温度、压力都发生突变。这一变化并不是均匀传播至空气中的每一部分的，因此某一时刻某一位置的压力变化不能代表整个药室内的压力变化。
图6为根据不同时刻的压力云图所选取的特征点，在同一压力范围内选择出一个最佳的特征点作为代表，一共选择了8个特征点,它们分别代表了不同时刻不同区域的压力变化情况。8个特征点的压力-时程曲线如图7所示。
图6 特征点所选位置
Fig.6 Selected positions of feature points
图7 8个特征点的压力-时程曲线
Fig.7 Pressure-time history curves of 8 characteristic points
由图7可知，曲线之间不存在重合的现象，因此所选的8个点没有多余点。 根据表1参数，计算85 mm口径某火炮的药室容积约为3 940 cm3，而测压器的体积约为22 cm3。根据国军标要求，测压器所占体积不可以超过药室容积的2.5%。依据本文火炮药室的参数计算，测压器所放个数不可以超过4个，而上述分析特征点共有8个，超过了测压器所能放置的个数，因此需要对这些点进行取舍。从曲线图可以看出，A、B、E、G这4个点分别在药室的底部、中间部分以及顶部。曲线出现峰值的时刻，以及峰值大小差异明显，因此选择这4个特征点来放置测压器。
3 数值分析
炸药在无限大的区域内爆炸，影响空气爆炸冲击波波面超压Δp的因素主要有4个因素，可以表示为：
Δp= f (E0,p0,ρ0,R)
(3)
式(3)中：E0为炸药的能量；p0为空气初始状态下的压力； ρ0为密度；R为传播的距离。
而炸药的能量E0可用其质量W与爆热Qν的乘积表示，因此可以将其简化为：
(4)
而常用的自由场爆炸冲击波超压峰值计算式有：
(5)
(6)
(7)
(8)
式(5)—(7)中： Δp为冲击波超压峰值，
为比例距离，
为测点爆心距，m；WT为TNT的装药质量，kg。
式(5)—(7)所预测的压力值是早年根据实验结果所得出的，而且是在自由场区域中的结果，本文的仿真实验是在密闭空间中进行，而且空间范围比较狭窄，因此根据仿真结果，重新进行了拟合，拟合公式为：
(9)
图8为上述公式为R在0.2～0.75 m·kg-1/3的压力变化曲线，由图8可以看出，式(9)跟仿真结果的契合度最高。
图8 各预测公式比例距离-压强曲线
Fig.8 Proportional distance-pressure comparison of each prediction formula
文献[15-16]中研究了自由场火药爆炸的压力时程曲线，可以由如下方程进行描述，即：
升压：
(10)
降压：
(11)
式(10)—(11)中:升压时间Tr=7.994 6×10-5R0.85；降压时间Td=7.167 3×10-4R1.342 5W0.313 6；控制衰减率的常数
由于本文仿真实验与文献[15]的仿真环境有所差异，故在升压时间与降压时间的取值上有所差异。它们二者的值是根据在药室内取多个点，分别观察不同点的升压时间与降压时间，取其均值，最终取Tr=30 μs，Td=150 μs。
图9为F点的局部压力曲线与式(10)—(11)的压力-时程曲线，将F点的压力曲线移至零点，2条曲线的拟合情况较好。其余的3个特征点的压力-时程曲线也可以用同样的方式获取。
图9 F点局部以及拟合公式压力-时程曲线
Fig.9 The partial pressure-time history curve of point F and the fitting formula
4 结论
本文中虽未考虑弹丸挤进过程，所得的压力曲线与真实膛压曲线有一定差别，但在弹丸挤进之前，曲线吻合度较高，选择不同位置处的点，观察膛压的变化情况，得出了以下结论：
1) 药室内不同位置处的压力数值变化存在很大差异，因此不可以用某一点位置处的压力变化情况近似为整个膛内的压力变化；
2) 依据不同时刻的压力云图，对于85 mm口径的某火炮而言，在考虑国军标要求以及经济性的前提下，共需放置4个测压器，而其余4个特征点可以用拟合的方式得出压力-时程曲线，以便完整的描述膛内压力变化情况。
[1] 郭天吉,吴宏斌.火炮内弹道环境测试技术[J].电子技术与软件工程,2019(21):115-116.
Guo T J,Wu H B.Ballistic environment test technology in artillery[J].ElectronicTechnology & Software Engineering,2019(21):115-116.
[2] 贺磊,姚养无,李树军,等.基于改进粒子群算法的火炮内弹道多参数符合计算[J].火力与指挥控制,2021,46(11):165-169.
He L,Yao Y W,Li S J, et al.Multi parameter coincidence calculation of gun interior ballistics based on improved particle swarm optimization algorithm[J].Fire Control & Command Control，2021，46(11)：165-19.
[3] Wang D Y,Nan F Q,Du P,et al.Simulation and calculation of interior ballistic of a large caliber artillery[C]//IOP Conference Series:Materials Science and Engineering,2019 International Conference on Optoelectronic Science and Materials 20-22 September 2019,Hefei,China,2019.
[4] Zhao J,Xu B,Wang Z C,et al.Modeling and simulation of interior ballistics based on actual combustion characteristics of propellants[J].Journal of Energetic Materials,2015,33(03):155-166.
[5] Zhang X B,Cheng C.Modeling of interior ballistic gas-solid flow using a coupled computational fluid dynamics-discrete element method[J].Journal of applied mechanics,2013,80(03):185-190.
[6] 裴东兴,沈静华,张瑜,等.放入式火炮膛压测试技术研究[J].火炮发射与控制学报,2017,38(03):69-73.
Pei D X,Shen J H,Zhang Y,et al.Research on bore pressure test technology of drop-in artillery[J].Journal of Gun Launch & Control,2017,38(03):69-73.
[7] 舒跃飞,李新娥,刘雪飞,等.瞬态高压多点压力测试系统[J].中国测试,2016,42(03):67-72.
Shu Y F,Li X E,Liu X F,et al.Transient high pressure multi-point pressure test system[J].China Measurement &Test,2016,42(03):67-72.
[8] 舒跃飞.基于光触发的微小型电容式测压系统研究[D].太原：中北大学,2016.
Shu Y F.Research on micro-miniature capacitive pressure measuring system based on light triggering[D].North University of China,2016.
[9] 王亚军.火炮膛内多点压力同步触发测试方法研究[D].太原：中北大学,2013.
Wang Y J.Research on multi-point pressure synchronous trigger test method in artillery chamber[D].North University of China,2013.
[10] 刘小龙,裴东兴,李新娥,等.火炮膛内多测试装置同步测试方法[J].中国测试,2017,43(12):130-134.
Liu X L,Pei D X,Li X E,et al.Simultaneous test method for multiple test devices in artillery chambers[J].China Measurement & Test,2017,43(12):130-134.
[11] 曾晓虹.爆炸冲击荷载下隔震结构连续倒塌模拟[D].兰州:兰州理工大学,2017.
Zeng X H.Simulation of continuous collapse of seismic isolation structures under explosive impact loads[D].Lanzhou University of Technology,2017.
[12] Sahu A,Omar Y A D,Sauer R A,et al.Arbitrary Lagrangian-Eulerian finite element method for curved and deforming surfaces[J].Journal of Computational Physics,2020,407:109253.
[13] Hwang Chan Gyu.Simulation of dynamic fracture and fluid-structure interaction in solid propellant rockets:Part 1 (theoretical aspects)[J].Journal of the Korea Academia-Industrial cooperation Society,2008,9(02):286-290.
[14] 尹建平,王志军.弹药学[M].北京：北京理工大学出版社，2012.
Yin J P,Wang Z J.Ammunition[M].Beijing:Beijing Institutionof Technology Press,2012.
[15] 杨亚东,李向东,王晓鸣.爆炸冲击波空中传播特征参量的优化拟合[J].爆破器材,2014,43(01):13-18.
Yang Y D,Li X D,Wang X O.Optimal fitting of characteristic parameters of airborne propagation of explosion shock waves[J].Explosive Materials,2014,43(01):13-18.
[16] 曹涛,孙浩,周游,等.近地爆炸冲击波传播特性数值模拟与应用[J].兵器装备工程学报,2020,41(12):187-191.
Cao T,Sun H,Zhou Y,et al.Numerical simulation and application of shock wave propagation characteristics of near-earth explosion[J].Journal of Ordnance Equipment Engineering,2020,41(12):187-191.