某大口径枪械内弹道挤进时期的数值仿真与分析

陆 野，周克栋，赫 雷，冯国铜，李峻松

（1.南京理工大学 机械工程学院，南京210094；2.中国兵器第二○八研究所，北京102202）

自动武器内弹道学是研究弹丸在膛内运动规律的科学，是枪械和弹药设计的理论基础，而弹丸挤进时期是影响内弹道初始条件的重要因素，对弹丸初速有重要影响，对研究弹丸与内膛的磨损关系、研制新型枪弹以及改进现有枪弹的射击精度具有重要的现实意义［1－2］，弹丸挤进经历的时间极短，难以用实验准确测量，因此经典内弹道学中通常假设在达到挤进压力瞬间弹丸才开始运动，忽略了弹丸挤进这一重要过程。近年来，国内外对此进行了多方面的研究［3－5］。本文以弹丸挤进过程为研究对象，采用多软件综合仿真的方法，模拟了弹丸的挤进过程，着重对枪管与弹丸挤进过程中的应力以及变形情况进行了数值仿真，获得了弹丸的嵌入阻力，通过回归分析获得了嵌入阻力的最优拟合公式。

1 弹丸挤进枪管的有限元建模

1.1 几何模型的建立

本文以某大口径机枪枪管及弹丸为研究对象，如图1所示，枪管壁厚、膛线等均未做简化，内膛分为弹膛、坡膛和线膛，由于仅分析弹丸挤进阶段，因此对枪管长度进行了缩短，弹丸由被甲、铅套、钢心等组成，弹丸内部以及枪管的构造如图1所示。

图1 弹丸挤进枪管三维模型

1.2 材料模型

枪管材料采用30SiMn2MoVA，弹丸被甲材料为覆铜钢片，弹丸内部由铅套包裹着钢心，由于弹丸和枪管材料在挤进坡膛的过程中出现高应变率、材料温度升高、应力软化等一系列复杂的非线性力学形态，本文采用了Johnson－Cook本构模型［6－7］。

Johnson－Cook本构模型可以写成如下形式：

式中：A，B分别为参考应变率和参考温度下材料的屈服强度和应变硬化系数；cn为参考应变率下材料的应变硬化指数；C为应变率敏感系数；cm为温度软化指数；εpl为变形 量为等效应变率为 参 考应变率；T为温度；Tr为室温；Tmelt为熔化温度。

1.3 网格精细划分

运用有限元专用前处理软件HYPERMESH对枪管和弹丸进行网格划分。如图2所示，枪管网格类型采用八节点八面体单元，其中对于膛线和坡膛等复杂结构处的网格进行了精细处理；弹丸采用SOLID164单元进行划分，其中对弹丸被甲的网格进行了精细处理。采用适用于大变形的单点缩减积分算法，避免了单元剪切闭锁问题；设置粘性沙漏控制，解决了沙漏问题，克服了零能模式；由于弹丸被甲单元会随着弹丸被甲的大变形而引起单元网格的大变形，甚至产生网格畸变，因此在计算过程中采用ALE（Arbitrary Lagrangian Eulerian）技术对网格进行更新，使整个计算过程中保持高的网格质量，确保数值计算的精度及效率。

图2 枪管膛线与弹丸网格划分示意图

1.4 弹丸挤进时期弹后压力的分析求解

由击发底火到弹丸全部嵌入膛线的过程称为挤进膛线时期，此过程由于弹丸移动量不大，可近似认为火药定容燃烧，火药气体定容燃烧状态方程为［8］

式中：f为火药力，Δ为装填密度，ψ为火药已燃烧的相对质量，ρm为火药密度，α为余容，pB为点火药压力。

燃气生成方程与燃烧速度方程为［8］

式中：Z为火药已燃烧的相对厚度；u1为燃速常数；e1为火药原厚度；χ、λ、μ均为药形系数；n为燃速指数。

联立方程（1）和方程（2），代入某大口径枪弹的火药特性参数，通过MATLAB编程即可计算出挤进膛线时期的弹底压力pd的曲线，由于自动武器挤进过程历时极短，本文进行仿真时仅取该曲线0～0.4ms区间，如图3所示，将此曲线施加于弹丸尾端面，即为力边界条件。设置枪管外表面的固定约束，即为位移边界条件。

图3 弹底压力曲线图

2 数值计算与仿真结果分析

本有限元模型采用显式非线性求解器ANSYS／LS－DYNA进行计算。针对弹丸被甲的变形与膛线刻痕的形成情况、弹丸挤进过程中应力、变形量以及嵌入阻力的变化情况，使用后处理软件Ls－Prepost，对仿真计算的数据进行整理分析，获得了弹丸挤进过程中的相关参数。

2.1 弹丸被甲膛线刻痕的形成过程分析

为更加直观地反应弹丸被甲膛线刻痕的形成过程，取膛线刻痕形成的阶段特征图，按时间顺序展示，如图4所示，随着弹丸挤进枪管，弹丸被甲膛线刻痕由弹头向弹尾延展，逐渐加深成形，且刻痕宽度逐渐变窄，不同膛线刻痕形成过程、形状基本一致。

图4 膛线刻痕的形成过程

2.2 弹丸所受应力、变形量分析

挤进过程中，弹丸与枪管之间的作用力主要有弹丸与膛线间的压力、摩擦力和导转侧力等，分析弹丸在挤进过程的应力、变形量的分布变化，对于枪管坡膛设计及弹丸设计具有重要的意义。

2.2.1 弹丸在X方向所受应力、变形量分析

如图5所示，为研究弹丸在X方向（由弹头指向弹尾）所受应力σ、变形量S的情况，取弹丸即将与膛线作用形成刻痕的单元（黑色单元区域）为分析对象，该组单元沿X轴方向分布且具有和膛线相同的缠度。

通过获取所取单元中沿X轴方向每一个单元的应力－时间曲线，绘制该组单元所受应力随X方向及时间变化的三维时空分布图，如图6所示，所绘制三维图选取仿真时间区间为0～0.4ms，测量弹丸刻痕沿X轴方向的长度，选取X轴坐标范围0～25mm。通过将计算得到的Mises应力云图三维数据化（见图6），更能直观表达其中的变化情况。

可以看出，弹丸在挤进瞬间应力就达到峰值183MPa，最大应力位置随着时间的增大，沿X轴方向向刻痕后端推移，由于弹丸结构设计为中间部分有3mm宽的凹造型，因此图6中X在12～14mm一段应力值较小；同时，在出现最大应力时刻后，应力值明显下降并趋于稳定。

图5 挤进前后X方向所选单元组示意图

图6 Mises应力随X方向及时间变化的三维时空分布图

以相同方法绘制所取单元的变形量随X方向及时间变化的三维时空分布图，如图7所示。

图7 弹丸被甲变形量随X方向及时间变化的三维时空分布图

可以看出，弹丸被甲在与枪管接触时刻变形量瞬间达到最大值1.18mm，然后变形量平稳无变化；弹丸被甲变形的发生时间随着弹丸的前进而增加，在时间－X轴平面呈斜线，整个挤进变形的时间区间为0.1～0.3ms；由于弹丸结构被设计为中间部分的凹造型，因此图7中X在12～14mm一段无变形。

2.2.2 弹丸圆周方向所受应力、变形量分析

为研究弹丸在圆周方向所受应力的时空分布图，取弹丸沿圆周方向的单元（即柱坐标θ＝0°～45°范围，图中黑色区域）为研究对象，如图8所示。

图8 挤进前后圆周方向所选单元组示意图

以圆周方向（θ）与时间为变量，绘制所选单元应力随圆周方向（θ）及时间变化的三维时空分布图，如图9所示。可以看出，在非膛线刻痕区，即θ≈0°～10°，35°～45°范围，弹丸所受应力值较小，与膛线刻痕区有明显区别；在膛线刻痕区，即θ≈10°～35°范围，单元应力在同一时刻出现峰值，由于枪管为右旋膛线，弹丸挤进时膛线刻痕左侧的应力值要大于右侧，图中显示为应力在膛线刻痕区逐渐变大50MPa，在θ≈30°处达到最大值187MPa。

图9 Mises应力随圆周方向（θ）及时间变化的三维时空分布图

以相同方法绘制所取单元的变形量随X方向及时间变化的三维时空分布图，如图10所示。可以看出，单元在非膛线刻痕区基本没有变形量，与膛线刻痕区有明显区别；所取单元在同一时刻瞬间达到最大值，而后趋于平稳；弹丸挤进时膛线刻痕左侧变形量比右侧大0.2mm，图中表现为变形量随圆周方向角度的增加而变大，在θ≈30°处达到最大值1.20mm。

图10 弹丸被甲变形量随圆周方向（θ）及时间变化的三维时空分布图

2.3 挤进过程中弹丸的嵌入阻力分析

弹丸挤进过程对于获得预期速度和射程具有重要的意义，挤进动作相当于延迟发射，从而增大了膛内平均压力，对于弹丸初速具有重要的影响。然而，目前自动武器行业对于挤进压力的定义较为模糊，在进行内弹道计算时，自动武器的挤进压力通常取为弹丸开始运动时的瞬间压力。目前国外对于挤进压力的定义主要有3种，其中应用最广泛的是将其定义为弹丸嵌入膛线的最大阻力。本文针对挤进过程中弹丸的挤进阻力，对数值仿真结果进行了拟合与比较，得出了可靠的弹丸嵌入阻力公式。

通过对仿真结果中每个单元所受的阻力进行耦合得到弹丸嵌入阻力，如图11所示。弹丸在0.1ms与枪管膛线接触，嵌入阻力瞬间达到峰值11.6kN；在0.3ms时完成挤进过程，嵌入阻力之后缓慢下降。

图11 弹丸嵌入阻力随时间变化的仿真计算结果

采用Gauss法、Sina法、Lorentz法、Voigt法等4种方法对弹丸嵌入阻力曲线进行拟合，如图12所示，通过计算每一条拟合曲线的可决系数R2，分析其拟合优度，得到最优的拟合曲线和公式［8］。

可决系数R2是测定多个变量间相互关系密切程度的统计分析指标，其值越大说明样本观测点与回归线靠得越近，模型拟合优度越好，其计算公式为

式中：为根据拟合方法所采用的回归方程计算出的估计值，¯F为样本平均值。

图12 4种方法得到的弹丸嵌入阻力拟合曲线

通过计算，Gauss法、Sina法、Lorentz法、Voigt法的可决系数R2分别为0.869 21，0.873 38，0.866 5，0.858 31。通过比较可知，最优的拟合方法为Sina法，得到拟合出的弹丸嵌入阻力公式为

式中：t为弹丸运动的时间，单位为 ms；Y0，tC，A，w为拟合常数。

3 结论

本文针对某大口径枪械挤进内弹道时期建立了弹丸挤进枪管的有限元模型，在此基础上进行了数值仿真，通过数据整理与分析计算，得出以下结论。

①精确、直观地模拟了弹丸被甲膛线刻痕的形成过程。随着弹丸挤进枪管，弹丸被甲膛线刻痕由弹头向弹尾延展，逐渐加深成形，且刻痕宽度逐渐变窄。

②弹丸所受应力的时空演化规律：在X方向上，刻痕区θ≈22.5°处接触应力在0.26ms达到峰值183MPa，最大应力位置随着时间、沿X轴方向向刻痕后端推移；在圆周方向上，由于枪管为右旋膛线，弹丸挤进时膛线刻痕左侧（导转侧）的应力值要大于右侧（非导转侧）50MPa，应力在膛线刻痕区逐渐变大，在θ≈30°处达到最大值187MPa。

③弹丸被甲的变形量时空演化规律：在X方向上，刻痕区θ≈22.5°处变形量在0.26ms达到最大值1.18mm，之后变形量平稳无变化；弹丸被甲变形的发生时间随着弹丸的前进而增加，在时间－X轴平面呈斜线，整个弹丸被甲挤进变形的时间区间为0.1～0.3ms；在圆周方向上，变形量随圆周方向角度的增加而变大，膛线刻痕左侧（导转侧）变形量比右侧（非导转侧）大0.2mm，在θ≈30°处变形量达到最大值1.20mm。

④针对数值仿真得到的弹丸嵌入阻力曲线，采用4种不同方法进行拟合，引入可决系数R2，得出了最优的嵌入阻力计算公式。

［1］张喜发，卢兴华.火炮烧蚀内弹道学［M］.北京：国防工业出版社，2001.ZHANG Xi－fa，LU Xing－hua.Erosion interior ballistics of gun［M］.Beijing：National Defense Industry Press，2001.（in Chinese）

［2］金志明.枪炮内弹道学［M］.北京：北京理工大学出版社，2004.JIN Zhi－ming.Interior ballistics of gun［M］.Beijing：Beijing Institute of Technology Press，2004.（in Chinese）

［3］樊黎霞，何湘玥.弹丸挤进过程的有限元模拟与分析［J］.兵工学报，2011，32（8）：963－969.FAN Li－xia，HE Xiang－yue.Finite element simulation and process analysis of projectile entering into barrel［J］.Acta Armamentarii，2011，32（8）：963－969.（in Chinese）

［4］韩文祥.弹丸挤进变形模拟问题研究［D］.南京：南京理工大学，2009.HAN Wen－xiang.The simulation research on the deformation process of projectile pulling into bore［D］.Nanjing：Nanjing University of Science and Technology，2009.（in Chinese）

［5］葛建立，杨国来，陈运生，等.基于弹塑性接触／碰撞模型的弹炮耦合问题研究［J］.弹道学报，2008，20（3）：103－106.GE Jian－li，YANG Guo－lai，CHEN Yun－sheng，et al.A study on projectile－barrel coupling problem based on elastoplastic contact／impact model［J］.Journal of ballistics，2008，20（3）：103－106.（in Chinese）

［6］JOHNSON G，COOK W.Fracture characteristics of three metals subjected to various strains，strain rates，temperatures and pressures［J］.Engineering Fracture Mechanics，1985，21：31－48.

［7］JOHNSON G R，COOK W H.A constitutive model and data for metals subjected to large strains，high strain rates and high temperatures［C］／／The 7th International Symposium on Ballistics.Netherlands：IBC，1983：78－96.

［8］李杰明.统计学原理［M］.上海：复旦大学出版社，2010.LI Jie－ming.General theory of statistics［M］.Shanghai：Fudan University Press，2010.（in Chinese）