尾翼弹高低空气动特性仿真及外弹道计算

弹箭气动参数对气动外形设计及弹道的计算具有重要意义，为此众多学者对弹箭气动参数的计算进行了大量的研究。文献[1]研究了大长径比超音速卷弧尾翼火箭弹在标准气象条件下,卷弧尾翼对全弹气动特性的影响；文献[2]通过研究得到某35 mm超音速榴弹在标准气象条件下阻力系数随弹形的变化规律；文献[3]研究了在标准气象条件下，某型导弹的亚音速气动特性并给出了该弹的气动参数；文献[4]研究了某小型单兵巡飞弹在2.5 km海拔高度处的气动参数，为巡飞弹的气动外形的设计提供参考；文献[5]研究了 4 km高空条件下，有无弹带对高速旋转弹丸气动特性的影响，得到弹带使得阻力系数有一定的升高，但弹带形状对弹丸升力系数影响不大的相关结论；文献[6]采用计算流体动力学对M910弹丸在全音速范围内静压为101 325 Pa,静温为 288 K条件下的气动特性进行了模拟，得到各气动参数与实验数据的误差在10%以内，具有较高的计算精度。文献[7]分析得到了某旋转弹丸在全音速范围内，静压为101 325 Pa,静温为292 K条件下，静态空气动力系数与大多数动力导数具有很好的一致性。文献[8]对来流马赫数Ma=0.7～2.7，攻角为 0°情况下的旋转弹丸绕流场进行了数值模拟,所得阻力系数结果与实验值和半经验公式计算值符合较好。文献[9]通过运用牛顿撞击理论计算得到某高超声速制导炮弹在静压为101 325 Pa，静温为293 K条件下的阻力系数与实验结果吻合较好。因此,该计算方法能够可靠地获得高超声速制导炮弹的气动力特性。文献[10]采用CFD方法，基于SST湍流模型对超声速火箭弹在特定海拔高度下的外流场进行了计算，计算结果表明该数值方法能较好地计算气动力参数。

然而在以往对弹箭气动参数的研究中，国内外很多学者只研究了固定气象条件下，弹箭的气动参数。少有人研究弹箭气动参数随不同海拔高度处的变化规律，但是由于不同海拔高度处的气象参数变化较大，则在不同海拔处发射的弹箭在飞行过程中受到的气动力也将会有较大的变化，从而影响弹箭的弹道。因此本文以小口径尾翼弹为研究对象，研究该尾翼弹随不同海拔高度，不同马赫数处的气动特性及其外弹道特性的变化规律，为小口径尾翼弹在不同海拔高度处的应用奠定基础。

1 数值计算方法

本文研究海拔高度在20 km以内，小口径尾翼弹的气动特性，其所属流域为连续流域，气体的间断效应并不显著，连续性假设成立，所以仍可运Navier-Stokes方程进行计算，选用雷诺平均的N-S方程作为控制方程。

控制方程为[11]：

式中：U为解向量；F和G为通量向量；Ma∞为来流马赫数；γ为比热比，Re∞为来流雷诺数。

湍流模型采用涡黏模型中的Shear-Stress Transport(SST) k-ω模型，求解2个方程的输运方程，即湍动能k,湍流频率ω，SST k-ω模型考虑到了湍流剪切应力的输运，能较好的模拟弹箭的气动力，因此选用该湍流模型计算尾翼弹的气动参数，k和ω的输运方程具体形式为[12]：

尾翼弹外形参数如图1所示。

头部母线形状为锥形；弹头部长径比

全弹长径比

前缘后掠角χ0=70°；尾翼平面形状为梯形；展弦比λW=LW/bav=1.2，Lw为尾翼展长，bav为尾翼平均几何弦长，bav=(br+bt)/2；根梢比

计算参考面积为弹体横截面积：S=7.85e-5m2。

综合考虑计算精度与计算效率，将计算域分层，在大的圆柱外流域内再增加一个小圆柱加密区域[13]，大圆柱形外流域直径为18倍弹径长度为10倍弹长，小圆柱流域直径为5倍弹径，长度为3倍弹长。采用求解黏性子层的方法对近壁面区域进行计算，即需要保证y+<5使得第一层网格节点位于粘性子层内[14]，并保证边界层数为8，划分得到如图2和图3所示的多面体网格。

由美国1976年标准大气分段逼近公式[15]计算：

计算位势高度H,它与几何高度Z的关系为：

式中：R0=6356.766 km,为地球平均半径；Z单位为 km。

1) 当0≤Z≤11.019 km时：

p/p0=W5.252 9, ρ/ρ0=W4.255 9

2) 当11.0191 km≤Z≤20.063 1 km时：

p/p0=0.119 53W，ρ/ρ0=0.158 98W

式中：p0=1 013.25 hPa，密度ρ0=1.225 kg/m3。

3) 黏性系数随高度变化的计算公式：

式中： βa=1.458×10-6kg/(s·m·K1/2)，Ts=110.4 K。

4) 理想气体定熵流动的当地声速计算公式为：

式中：k为绝热指数，空气k=1.4; Rg为气体常数，空气Rg=287 J/(kg·K)；T为热力学温度(K)。

由以上公式计算得到各海拔高度的气象参数值，作出相关曲线如图4所示。

设置流场外壁面为压力远场边界，弹体表面采用无滑移壁面边界条件，求解方法采用压力-速度耦合(Coupled)算法，对流通量采用二阶迎风格式进行离散。对气动参数计算时，远场的温度、压力、密度、黏性系数等气象参数已由以上部分计算得到。假设来流为理想气体，来流马赫数为Ma=0.6～3，攻角为α=0°。

算法有效性验证：

本文选用标准模型(Non-Rolling Basic Finner)进行数值算法有效性验证，该标准模型已由美军弹道研究实验室进行了大量风洞试验。图5为其几何外形参数，弹体直径d=20 mm，以弹体横截面积为参考面积,并用以上数值方法模拟标模弹的外流场，得到零升阻力系数，然后作出图6。

图6所示的仿真零升阻力系数曲线和文献[16]的实验值拟合曲线大致吻合，误差小于10%，在允许范围以内。由此可知，本文中所用的数值计算方法具有较高的可信度。

2 计算结果与分析

由数值方法计算得到不同海拔高度处的尾翼弹阻力值，作出图7：在同一马赫数下，尾翼弹阻力随着海拔高度的增加在减小，同一海拔高度处，尾翼弹阻力随着马赫数的增加近似线性增加。

由公式

计算得到在不同海拔高度处及不同马赫数下的阻力系数Cd值，并作出图8、图9，式中：D为尾翼弹阻力，ρ∞为来流密度，V∞为来流速度，S为弹体横截面积。从图8中可以看出阻力系数随着马赫数的增加先急剧增加，而后缓慢减小，并且在Ma=1.1附近取得最大值。这是由于弹箭在亚声速段，弹箭阻力由摩阻和涡阻组成，在跨声速飞行时，弹体表面产生弹体激波，弹箭阻力由摩阻、涡阻和波阻组成，在超声速段飞行时，弹头部脱体激波转变为附体激波，且由图10不同海拔高度不同马赫在α=0°时的速度云图可以看出，激波锥角随着马赫数的增大而减小，因而波阻减小，弹体阻力系数减小。

在Ma=0.6到Ma=3全声速范围内，阻力系数随着海拔高度的升高而增加，在同一海拔高度处，阻力系数随着马赫数的变化趋势相同。

由数值方法计算得到弹体各部分阻力系数变化值，作出图11所示曲线，并可以得到：随着海拔高度的增加，尾翼弹头部，弹体部，尾翼部的阻力系数都随着海拔高度的增加而增加，其中尾翼部的阻力系数对海拔高度最为敏感，而弹底阻力系数随海拔高度增加而降低。因此尾翼部的阻力系数升高是造成阻力系数随着高度增加的主要原因。

从图8曲线得到：在海拔高度17 km内，阻力系数与海拔高度基本为线性关系；亚音速范围内的阻力系数增加速率要高于超音速范围的增加速率；在海拔高度超过17 km后，亚音速段，跨音速段的阻力系数增加较快，而超音速段的增长速率基本不变，亚音速段，跨音速段的阻力系数对海拔高度的变化更为灵敏。

3 阻力系数拟合

文献[17]提出：用Logistic曲线与三次抛物线分俩段拟合亚音速段(1.1Ma为界)，超音速段弹丸阻力系数，该方法拟合误差较低且分段数较少，方便于外弹道编程使用。

拟合亚音速段Logistic曲线解析表达式为：

式中：a分别为Logistic曲线下渐进线纵坐标值；b为上下两条平行的渐进线之间的距离；c，d为拟合曲线待定系数；x为弹丸飞行马赫数；y为弹丸零升阻力系数。

拟合超音速段三次抛物曲线解析表达式为：

式中：x为弹丸飞行马赫数；y为弹丸零升阻力系数；A、B、C、D为拟合曲线待定系数。

用确定系数R-square来表征拟合的程度：

式中：

拟合数据和原始数据对应点的误差平方和；

原始数据和均值之间的差

由上式可知R-square越接近1，表明所拟合方程的拟合程度高。用Matlab拟合方程R-square值均为0.99以上，拟合各海拔高度处阻力系数方程如表1，拟合曲线如图12。

4 外弹道计算

射角θ0=10°；初速v0=900 m/s；t=0,x=0；y=0 km,10 km,17 km; vx=vvcosθ0；vy=v0sinθ0

由质点外弹道计算得到各海拔高度处发射该尾翼弹的外弹道特性曲线，质点外弹道方程组[18]：

由质点外弹道方程计算得到弹道高，弹丸存速，弹道倾角随海拔高度的变化值，作出图13。

由图13可以得到：在不同海拔高度处，同一射角，初速下发射的弹丸弹道轨迹，弹丸速度，弹道倾角有较大的变化。因阻力值随着海拔高度的增加而降低，所以弹丸存速随着海拔高度的升高而增加，弹道高随着海拔高度的增加而下降，弹道倾角随海拔高度的增加而下降。

5 结论

在全声速范围内，同一马赫处，阻力系数随着海拔高度的增加而增大，且基本为线性关系；尾翼部的阻力系数随着海拔高度升高是造成整弹阻力系数升高的主要原因；在同一海拔高度处，阻力系数随着马赫数的变化规律相同。

全声速段各海拔高度处阻力系数的拟合方程具有较高的拟合精度。

在相同射角，相同初速，不同海拔处发射的尾翼弹的弹道轨迹，速度，弹道倾角不同，因其阻力值随着海拔高度的增加而降低，弹丸存速随着海拔高度的增加而增加，弹道高，弹道倾角随着海拔高度的增加而下降，需考虑海拔高度对弹道参数的影响。

[1] 谢志敏，杨树兴，陈伟.大长径比卷弧尾翼火箭弹气动特性数值研究[J].固体火箭技术，2009,32(06)：596-599.

[2] 张熠,周克栋,郝雷.某35 mm超音速榴弹弹丸阻力系数分析[J].兵器装备工程学报，2017(11):59-64,72.

[3] 陈卫东，唐小平，曾奎，等.基于工程和数值方法的导弹气动特性计算[J].航空计算技术，2012,42(03)：1-5.

[4] 陶福兴，张恒，李杰.一种小型单兵巡飞弹的气动外形设计[J].弹箭与制导学报，2015,35(06):111-114,118.

[5] 孟鹏，陈红彬，钱林方，等.弹带对高速旋转弹丸气动特性影响的数值模拟[J].兵工学报，2017,38(12)：2363-2372.

[6] DESPIRITO J,HEAVEY K R.CFD Computation of Magnus Moment and Roll Damping Moment of a Spinning Projectile[C].AIAA-2004-2713,2004.

[7] SIDRA I.Silton,Navier-Stokes Predictions of Aerodynamic Coefficients and Dynamic Derivatives of a 0.50-cal Projectile[C].AIAA -2011-2030,2011.

[8] SILTON S I.Navier-stokes computations for a spinning projectile from subsonic to supersonic speeds[J].Journal of Spacecraft and Rockets,2005,42(2):223-231.

[9] 朱胤,王旭刚.多尾翼式高超声速制导炮弹气动布局设计与特性分析[J].兵器装备工程学报,2020,41(01):48-52.

[10] 陈东阳，ABBAS LAITH K,等.超声速旋转火箭弹气动特性仿真和分析[J].计算机辅助工程，2013,22(06)：64-68.

[11] 李小林，杨帆，傅建明，等.MICA导弹气动性能分析研究[J].空天防御,2019,2(02):9-15.

[12] 陈东阳.超音速旋转弹箭气动特性及流固耦合计算分析[D].南京：南京理工大学,2014.

[13] 戴凯.基于FLUENT的飞行器气动特性仿真研究[D].西安：西安工业大学，2015.

[14] 胡坤，胡婷婷，马海峰,等.ANSYS CFD入门指南：计算流体力学基础及应用[M].北京：机械工业出社,2018.

[15] 韩子鹏.弹箭外弹道学[M].北京:北京理工大学出版社，2014.

[16] MACALLISTER L C.The Aerodynamic Properties of a Simple Non Rolling Finned Cone-Cylinder Configuration Between Mach Number 1.0 and 2.5,Blallistic Research Laboratories Report No.934,1955.

[17] 杨翔，王雨时，闻泉.应用阻力系数拟合曲线解析式数值解算外弹道诸元[J].弹箭与制导学报，2014，34(05)：151-155.