受加工工艺、设备条件等因素的限制，实际生产的弹丸与理论设计会存在一定偏差。加工偏差的大小对高速飞行弹丸的气动特性与稳定性来说还是会产生一定的影响[1-4]。过大的偏差会直接影响弹丸散布和飞行稳定性，所以对加工误差造成影响的问题进行研究是十分必要的。陈东阳等[3]建立了考虑结构误差(包括质量偏心、弹体不同轴)的气动特性计算模型，通过与无结构误差模型的计算结果进行对比得到了受各结构误差项影响显著的气动力参数。S.I.Silton[5]利用CFD的方法对旋转弹丸从亚音速到超音速的气动特性参数进行了模拟计算。李建中等[6]分析了某型火炮弹丸的质量偏心对弹丸起始扰动的影响。Albert E等[7] 研究了弹丸的质量不对称对弹形系数的影响，其研究表明在弹丸基本形态都相同的情况下，主轴的偏差和重心偏移对弹丸的弹形系数造成的影响不容忽视。
由于数值模拟技术不仅可以模拟出极限条件，还能够反映弹丸飞行的基本规律，获取精细化的数据，因此成为研究弹丸飞行问题的重要手段之一，从而被广泛应用。本文选取典型的榴弹进行研究，设定弹径为82 mm，基于此采用数值模拟的方法进行仿真分析。在弹丸加工诸多误差中，质量偏心误差是不可忽视的因素之一，本文主要研究质量偏心误差单因素对飞行弹道产生的影响。针对无结构误差弹丸与质心偏差弹丸，采用SST k-ω湍流模型和滑移网格技术[8]处理弹头的高速旋转所引起的动边界问题，模拟了两种弹丸的气动特性。结合弹丸的气动参数和膛口扰动状态，利用弹丸6自由度外弹道方程对其外弹道过程进行求解，得到两种弹丸的攻角衰减曲线。并通过解析算法分别求解两弹丸的急螺稳定因子，动态平衡角和动态稳定因子的大小从而对比两种弹丸的飞行稳定性。
1 数值模拟
以三维Navier-Stokes方程(下文简称N-S方程)为基础，采用剪应力输运湍流模型(Shear Stress Transport k-ω,下文简称SST k-ω)分别对无结构偏差与质心偏差的旋转榴弹气动特性及外部流场分布情况进行数值模拟。利用solidwork建立3D物理模型并用ICEMCFD画出计算域网格。由于旋转榴弹出膛后具有较高的转速，为准确模拟弹头旋转导致的非定常流动对其气动特性的影响，本文采用滑移网格方法处理弹头高速旋转所引起的动边界问题。
滑移网格技术[9]要求存在一个外部固定区和包围弹体的内部运动区，2个区域之间具有一对交界面。交界面上的网格节点不需要重合，只需要在滑移交界面上进行数据插值，即可保证2个区域之间的通量守恒。且内部运动区的网格单元在运动过程不发生变形，因而滑移网格技术还有占用计算机内存少、计算精度高的优点。
1.1 湍流模型
采用剪应力输运SST k-ω湍流模型，该模型是标准k-ω模型修正后的两方程模型，其不仅能够较好地模拟边界层内的低雷诺数流动，还能模拟边界层外充分发展的湍流流动。并且能适应压力梯度变化的各种复杂物理现象。该模型对湍流粘性系数进行了修正，且考虑了湍流剪切应力从而不会对涡流黏度进行过度预测，因而能更好地适用于跨声速及超声速来流下的弹头扰流场的模拟。
SST k-ω湍流模型的数学表达式为
(1)
(2)
式中：k和ω分别为湍动能和耗散率；σk和σω分别为k和ω的湍流普朗特常数；
为i方向的平均速度；xi和xj分别为i和j方向的位置坐标； μ和μt分别对应空气的动力黏度和湍流黏度；Gk和Gω分别为由平均速度梯度和浮力影响引起的湍动能产生项； β′和β″为模型常数；Dω为正交发散项。
1.2 弹体模型及边界条件
模型采用直径为82 mm的无结构偏差的旋转榴弹及与原质心径向偏离0.2 mm的结构偏差榴弹。如图1与图2所示。图3为计算域结构示意图，计算域将分为3部分，从内到外分别为边界层、内部旋转区和外部固定区。经过多次试算[10]，内部滑移区边界选择弹尾距离边界为2D，其他部位贴合壁面距离为1D；外部固定区域，前段距离弹尖7D，圆柱部位直径为17D。
图1 无结构偏差弹丸3D物理模型示意图
图2 质心偏差为0.2 mm弹丸模型示意图
图3 计算区域结构示意图
全弹面网格及全计算域空间网格分别采用正交性较好的四边形网格和六面体网格，网格数量约为34万。
采用的计算条件为：外部固定区的外边界采用压力远场边界条件(Pressure Far-Field)，总温T=300 K,总压P=101 325 kPa。弹丸攻角α为4°。内部旋转区与弹头壁面相关联，随壁面以相同角速度旋转；内部旋转区与外部固定区通过交界面传递数据，交界面采用滑移边界条件。选取压力远场进行初始化。如图4与图5为网格划分示意图。
图4 旋转区域网格示意图
图5 压力远场外边界网格示意图
选用的湍流模型是剪应力输运SST k-ω湍流模型。由于本文中的流体具有可压缩性，因此选取密度基求解器(Density-Based)。弹丸处于旋转状态，因此选用瞬态计算，并采用显式算法以节省计算机内存；差分格式选用高精度的二阶迎风格式。由于气动参数以系数方式输出，因此，需要先选定合适的参考值[11]。本文的数值模型需要设定的参考值有参考面和参考长度；其中参考面选取弹丸的横截面面积0.005 28 m2，参考长度为弹丸直径82 mm。
1.3 计算结果
在马赫数从0.3到0.8，攻角为4°的工况下通过数值模拟分别得出两弹丸的阻力系数、升力系数、俯仰力矩系数、滚转阻尼力矩系数和马格努斯力矩系数，计算结果如图6～图10所示。
图6 阻力系数曲线
图7 升力系数曲线
图8 俯仰力矩系数
图9 滚转阻尼力矩系数导数曲线
图10 马格努斯力矩系数导数曲线
从图6可知，无结构误差与质心偏移的弹丸阻力系数几乎相同，表明阻力系数的大小只和外形有关，与质心偏移无关。图7和图8分别是对升力系数和俯仰力矩系数的模拟结果，由图可知在弹丸低速飞行时，质心偏移对升力的影响不大。由图9和图10可看出，质心偏移对作用在弹丸上的滚转阻尼力矩系数和马格努斯力矩系数都会产生较大的影响，特别是对马格努斯力矩系数的影响尤为明显。
2 质量偏心对弹丸攻角衰减的影响
为了保证弹丸的稳定飞行，其飞行的攻角应当满足快速衰减的要求[12]，并最终保持在一个很小的波动范围之内。将本文模拟得到的气动参数代入六自由度外弹道程序可模拟出两弹丸在飞行过程中攻角的衰减过程，如图11、图12所示，由此可直观反应质量偏心对弹丸飞行稳定性的影响。
图11 无结构误差弹丸攻角曲线
图12 结构误差弹丸攻角曲线
由图11与图12两弹丸的攻角-时间曲线图可以看出，由于弹丸质量偏心的存在，弹丸攻角的衰减速度变慢，尽管其最终还是呈稳定状态，但是难以满足快速衰减的要求，势必影响弹丸的外弹道性能，对空间散布产生不利影响。因此弹丸的质量偏心必须控制在一定的范围内。
3 稳定性计算
3.1 弹丸急螺稳定性计算
弹丸相关参数如表1所示。
表1 弹丸相关参数
类别无结构偏差弹丸质心偏差弹丸极转动惯量Jx/(kg·m-2)0.002 942 3540.002 942 431赤道转动惯量Jy/(kg·m-2)0.008 8360.008 836弹丸全长l/mm179.3179.3弹丸头部长度lt0/mm4545质心到弹丸头部界面的距离h0/mm106106弹丸速度v/(m·s-1)238238弹径d/mm8282弹重m/kg2.32.3
弹丸飞行时，为保证稳定飞行[13]，其章动角δ应维持在有限的范围内变化。 在理论上，只要急螺稳定系数大于零，即可使弹丸保持稳定。由于射击精度的要求，就需要限制弹丸章动角的最大值。因此需要控制急螺稳定系数(σ)或者急螺稳定因子(S)的大小，其具体值一般通过试验确定。对于现代火炮弹丸系统中的榴弹[14]，急螺稳定系数取值一般为0.3～0.6。如果σ值取值过高，则将会使追随稳定性变差。因此σ的取值同时还要考虑弹道曲线段上的追随稳定性。
弹丸在直线段的飞行稳定性是通过其急螺稳定性来衡量的。本文通过求解急螺稳定因子S来判断其是否满足稳定性条件，当S>1时即满足条件。通过表1已知弹丸的参数。根据外弹道学[15]，弹丸的急螺稳定因子可以表示为
(3)
式中： η为火炮缠度，本例可按下式确定
a=πd=3.14d
b=a/tan4°≈44.8d
η=b/d=44.8
h/d可按下式计算：
h/d=h0/d+0.57(lt0/d)-0.16=1.445
kmz(M)的值取决于弹丸的全长l和初速v0，它可按以下公式计算
(4)
计算弹丸的急螺稳定性时，由于在炮口处弹丸的急螺稳定性最差，因此，通过求解炮口处的急螺稳定因子进行判断。两弹丸的计算结果如表2所示。
表2 计算结果
类别无结构误差弹丸质量偏心弹丸极转动惯量Jx/(kg·m-2)0.002 942 3540.002 942 431赤道转动惯量Jy/(kg·m-2)0.008 8360.008 836急螺稳定因子S2.143 62.143 7
由以上计算可看出，两弹丸的急螺稳定系数大小几乎是相同的，并且两弹丸都满足急螺稳定性条件。这是由于0.2 mm的质心偏差引起转动惯量的改变甚微，使得上述计算得到的结果差别很小。这表明质心偏差为0.2 mm的结构误差对弹丸起始运动时的急螺稳定性影响不大。
3.2 追随稳定性计算
弹丸在曲线段飞行时，弹轴的进动运动较之于直线段上的进动运动具有某些质的差别。此时，弹轴不再是绕弹道切线作圆锥运动，而是绕某一动力平衡轴作圆锥运动。动力平衡轴偏离弹道切线的夹角，称为动力平衡角，以δp表示。δp为弹丸的追随稳定性特征数。δp值越小，弹丸的飞行定向性也越好。δp在弹道上是变化的。但在弹道顶点附近，因速度翻转力矩的速度函数值等均达到最小值，而cosθ值最大，故相应的δp也最大。亦即该处的追随稳定性最差，所以应把此处的δp控制在符合要求的范围内，即小于某一允许值。顶点的飞行动力平衡角为
(5)
本文计算弹丸飞行射程为400 m时的运动情况。
式中：vs为弹道顶点的速度，根据外弹道程序得到
为弹道顶点的高度，
为取决于弹道高的函数，
综上
δps=0.104°
计算结果显示两弹丸在弹道顶点处都具有追随稳定性，且大小也相同，表明0.2 mm的质心偏差对弹丸的追随稳定性几乎没有影响。
3.3 动态稳定性计算
在考虑弹丸的动态稳定性时，不仅要考虑翻转力矩的作用，而应该综合考虑在所有力矩(包括翻转力矩，马格努斯力矩，赤道阻尼力矩)的作用下弹丸的攻角在弹道上的衰减情况。
为使得弹丸满足动态稳定，急螺稳定因子S和动态稳定因子Sd间必须满足下述条件，即
(6)
动态稳定因子可用下式来表示：
(7)
(8)
式中：
为升力系数导数；
为马格努斯力矩系数导数；
为赤道阻尼力矩系数导。
计算结果如表3所示。
表3 计算结果
类别质量偏心弹丸无结构误差弹丸升力系数导数C'y1.750 61.799 6马格努斯力矩系数导数m'y-0.043 5-0.083 2赤道阻尼力矩系数导数m'zz-0.000 342-0.000 419动态稳定因子Sd0.7380.515
从两弹丸的计算结果可看出，结构误差弹丸的动态稳定因子大于无结构误差弹丸的动态稳定因子，两弹丸都满足上述动态稳定性条件。
从以上3种稳定性的计算结果可看出0.2 mm的质心偏移虽然对弹丸的起始处和最高点处的稳定性影响微弱，但对弹丸的动态稳定性有较大的影响。这是由于质心偏移造成旋转轴线的偏移，从而导致气动特性产生很大的变化，最终对弹丸动态稳定性产生较大的影响。
4 结论
1) 旋转榴弹质心偏移的影响主要表现为使马格努斯力矩和滚转阻尼力矩发生变化，从而对弹丸的飞行稳定性产生影响。
2) 质量偏心会导致弹丸攻角的衰减速度变慢，增大弹丸的空间散布。弹丸质心偏移必须控制在一定的范围内。
3) 0.2 mm的质心偏移对弹丸的急螺稳定性和追随稳定性几乎无影响，而对动态稳定性造成较大影响。
4) 本文方法可以对质心偏移的影响程度进行定量分析计算，实现对弹丸的辅助设计和加工误差限的确定。
[1] 苗瑞生，吴甲生.旋转弹丸空气动力学[J].力学进展,1987,17(4):479-486.
[2] 赵子华.不对称因素引起射弹散布的分析与估算[J].弹道学报，1993,4(18):45-49.
[3] 郭锡福.质心位置对旋转稳定弹丸稳定性的影响[J].南京：华东工程学院学报，1982,1(007):81-95.
[4] COOPER G.Flight Stability of an Asymmetric Projectile with Activting Canards[J].Journal of Spacecraft and Rockets,2012,49(1):130-135.
[5] 陈东阳,ABBAS L K,芮筱亭,等.结构误差对旋转稳定弹丸气动特性影响的数值模拟[J].空气动力学报,2014,32(5):705-711.
[6] SILTON S I.Navier-Stokes Computations for a Spinning Projectile from Subsonic to Supersonic Speeds[J].Journal of Spacecraft and Rockets,2005,42(2).
[7] 李建中,胡敬坤,等.某型火炮弹丸质量偏心对弹丸起始扰动的影响分析[J].兵工自动化,2017,36(11):8-11.
[8] ALBERT E,HODAPP.Effect of Mass Asymmetry on Ballistic Match of Projectiles[J].Journal of Spacecraft and Rockets,1976,13(12):757-760.
[9] 纪兵兵,等.ANSYS ICEM CFD基础教程与实例讲解[M].北京：机械工业出版社,2015.
[10] 马杰.超声速旋转弹箭气动特性与流体分离控制[D].南京:南京理工大学,2016.
[11] 于勇,等.FLUENT入门与进阶教程[M].北京：北京理工大学出版社,2008.
[12] 王时雨.旋转弹丸外弹道自转角速度衰减规律半经验公式[J].探测与控制学报,2003，25(1):1-6.
[13] 吕铁刚，张亚,等.低转速条件下的弹丸飞行稳定性分析[J].现代防御技术,2018,46(3):61-65.
[14] 魏惠之等.弹丸设计理论[M].北京：国防工业出版社,1985.
[15] 韩子鹏.弹箭外弹道学[M].北京:北京理工大学出版社,2008:23-60.