低速旋转导弹滚转特性数值模拟研究

徐科杰，王学德，张 超，张隽研，王锌晨

(南京理工大学 能源与动力工程学院， 南京 210094)

现代制导武器中滚转控制是非常重要的问题，无论是带有控制系统的战术导弹、火箭弹、炮弹都需要控制其滚转速度，以此测算出飞行滚转角达到导弹飞行姿态控制的目的，从而提高导弹飞行稳定性和打击精度。导弹的稳定转速(即平衡转速)与导转力矩以及滚转阻尼力矩息息相关，因此研究旋转导弹的滚转特性极其重要[1]。

从20世纪70年代开始，国外学者利用风洞试验对导弹滚转特性进行了大量研究。F.J.Regan、Leroy、Collins.J.A[2-4]等对标准模型进行了风洞实验并得到了一些有用结论。随着计算机性能的大幅度提高，研究人员开始采用数值模拟方法对滚转特性进行研究。Pual Weinacht、Ameer G.Mikhail等[2-4]验证了基于N-S方程的数值计算方法计算导弹气动特性。Erdal Oktay、John F.Stalnaker等[6,7]采用基于不同算法的求解器对标准模型Basic Finner进行了非定常计算。在国内，邓帆[8]发现在超声速条件下的气流雍塞是导致栅格翼导弹的滚转阻尼导数减小的主要原因。周德娟[1]分别采用风洞试验和数值仿真对不同翼身组合导弹滚转特性进行了系统研究，两者计算结果吻合较好。余奇华、敬代勇[9]对不同转速下的旋转尾翼鸭式布局弹箭的流场特性进行了数值模拟，发现尾翼旋转可以显著提高其滚转控制能力。

国内外学者从风洞试验与数值仿真两方面对导弹的滚转特性进行了大量研究，基本上以研究马赫速、攻角与导弹滚转特性之间的关系为主，但在滚转特性随转速以及尾翼斜置角的变化方面研究较少。本文基于滑移网格技术，利用计算流体力学(CFD)软件对不同尾翼斜置角以及转速下导弹三维绕流流场进行了数值模拟，从流场结构以及气动特性两部分对导弹的滚转特性进行了阐述。通过数值模拟计算不同转速下的滚转力矩系数，使用Matlab进行数据拟合得到了滚转力矩系数随转速变化的一次拟合函数，进而得到导弹在不同工况下的平衡转速，并且与工程算法得到的平衡转速进行对比，得到快速计算导弹的滚转力矩系数以及平衡转速的经验公式。

1 数值模拟方法

滑移网格技术是20世纪80年代逐渐兴起的一种基于非稳态思想的数值模拟方法，近年来得到了广泛的应用，尤其是在弹箭旋转的气动仿真领域。图1为滑移网格技术的信息交换示意图。整个流动区域被分为包含1、2、3单元和4、5、6单元的两个不同的区域，其中一个区域相对另一个区域做相对滑移运动来模拟真实的流动过程。ABCD和EFGH是各自子区域与滑移面重合的边界面。在数值模拟过程中，1、2、3单元和4、5、6单元在运动过程中，会形成abcdefg节点，并沿着交界面做相对滑移运动。其中子块2的信息通过面cd和面de传递到子块4和子块5，其他单元之间的信息传递与此类似。

滑移网格技术作为一种特殊的动网格，适用于周期性旋转运动物体的流场计算，具有计算精度高，占用内存少等优势[10]。当模拟弹箭旋转运动时，只需要一个非旋转区域和一个包围弹体的旋转区域，两者之间存在一个交界面，交界面的网格节点数不必相同。仅仅在滑移网格交界面进行数据插值，即可保证两个区域之间通量守恒。图2(a)和图2(b)给出了网格区域划分示意图。

采用基于滑移网格下的三维积分守恒型非定常N-S方程组[11]，由质量守恒方程、动量守恒方程和能量守恒方程组成。

式中：ρ,U,e分别表示气流的密度、速度和单位体积的能量；vm为弹体的运动速度；P为静压力；I为单位张量；τ为黏性张量；q为热通量；S和v分别表示控制体的体积分和面积分的积分区域；r为S的外法向单位向量。

采用有限体积法对空间进行离散，对流项采用二阶通量差分分裂(FDS) AUSM+格式，黏性项采用中心差分格式。

Realizablek-ε模型是对RNGk-ε湍流模型做了修正而产生的，因此可以更为准确地反映流体真实的流动过程，尤其是旋转流动计算中气流的流动状态。

2 算法有效性验证

2.1 验证模型

选取美国陆-海军动导数计算标准模型(ANF)进行数值算法的有效性验证，ANF模型经过大量风洞实验和飞行试验验证，被广泛采用为动导数的计算标准模型，如图3所示，ANF-1与ANF-2模型的差别是：ANF-1无尾翼斜置角而ANF-2有尾翼斜置角。表1和表2给出计算条件，其中攻角(α)、马赫数(Ma)、雷诺数(Re)、总压(P0)、静压(P)和总温(T0)、静温(T)同风洞实验条件,Ω·为对应的导弹无量纲旋转速度(Ω·的定义为Ω·=D*Ω/V,D为弹体直径,V为当地声速，Ω为导弹转速(rad/s),N为外循环中每转动一周的总计算的迭代次数,Δt为外循环的计算步长,i为每一次外循环计算中内循环的迭代次数。

MaRe(106)P0/PaT0(k)Ω·2.491.8644 816311.10.034 6NΔt(10-6)iD/mm1 4401.1312031.75

表2 ANF-2模型的计算条件

利用ICEM商业软件，采用多块对接网格生成技术对整个计算区域生成六面体网格。纵向采用C型网格，横向采用O型网格。第一层距离物面为1.0×10-5m，网格伸展比为1.1，控制第一层网格到物面的距离和网格伸展比的目的是为了保证y+<1。图4(a)和图4(b)分别是ANF-1模型弹体头部和尾部附近网格。

弹箭表面设置为无滑移壁面条件，弹体跟随旋转区域以相同的转速同步运动。非旋转区域和旋转区域的交界面设置为滑移边界条件；非旋转区域的外边界设置为压力远场条件。

2.2 算法验证分析

图5(a)和图5(b)分别给出了不同攻角时ANF-1模型的滚转力矩动导数和俯仰力矩系数的计算值与实验值[12]。从图5可以看出滚转力矩动导数和俯仰力矩系数与相应的实验值吻合较好，最大误差分别为17.6%和8.62%(最大误差均出现在攻角α=80°附近)。同时发现攻角越大，误差越大。这是由于攻角越大，自旋导弹附近的流场变得更为复杂。为了进一步验证该数值模拟方法的精确性，对比了ANF-2模型的轴向力系数和法向力系数随马赫数的变化规律，并且与实验值进行了比较，如图6所示轴向力系数与法向力系数与实验值[13]吻合较好，其最大误差分别为3.3%(1.1马赫数下)和7.2%(3.5马赫数下)。可知该数值模拟方法在求解小攻角、全马赫数情况下有较高的精确度，验证了算法的正确性。

3 某低速尾翼弹滚转特性计算

计算模型为某尾翼弹，由锥形部、圆柱段和尾翼组成，其尺寸外形如图7所示,弹径为0.4 m，弹长为L为7.3D。选取表3中的参数作为计算状态，气动力和力矩计算的参考点位于弹体质心位置，参考面积为弹体最大横截面积，参考长度为弹径D。

3.1 网格尺度以及时间步长独立性验证

模型网格数目以及物理时间步长的选择对仿真计算结果有较大影响，尤其是在非定常计算中。若网格数目较少，将无法准确地捕获到流场信息甚至导致计算出错；若网格数目过大，则消耗计算资源，延长计算周期。与网格独立性原理一样，如物理时间步长太大，则会导致计算发散，如物理时间步长过小的话，又会加大计算时间。因此需要选取一个合理的时间步长以及网格数目。表4给出了Ma=2.0、α=2°、Ω·=0.015 1、γ=0°下，不同网格数下阻力系数CD、升力系数CL、俯仰力矩系数Cm以及滚转力矩系数Cl参数的计算结果。可见随着网格的不断加密，该导弹的气动特性参数逐渐趋于收敛。考虑到900万与800万网格下的计算结果差异较小(以900万网格下的气动数据作为基准，四组气动力参数的最小差值分别为1.3%、1%、1.2%、3.8%)，因此后续计算以800万网格为准。图8给出了导弹四组气动特性参数随物理时间步长Δt的变化曲线(计算条件参考表3所列的计算参数) ，考虑到Δt=0.000 4 和Δt=0.000 6中各气动特性参数差异较小(以Δt=0.000 4下的计算结果作为基准，四组气动力参数的最小差值分别为0.5%、1.1%、1.6%、1.3%)以及需要消耗的计算周期，本文采用Δt=0.000 6。

表3 模型的计算参数

表4 网格数对气动特性参数的影响

3.2 流场分析

如图9(a)和图9(b)所示为尾部截面的速度矢量图，从图9(a)可知，由于攻角的存在，使得在流场区域存在一个整体向上的速度分量。从图9(b)可以明显看出，导弹在旋转过程中，尾翼引起其附近气流随导弹一起顺时针旋转。A区域与C区域气流速度矢量方向相反，因此形成一个速度矢量为0的区域即B区域，同时也会在翼根位置即D区域形成一个较小的涡流。

3.3 气动特性参数分析

图10为马赫数Ma=2.0、尾翼斜置角γ=1°和不同转速下的尾翼位置的压力云图。从图10(a)可知，对于上下尾翼，由于整个尾翼处于背风(或者迎风)气流状态，当存在较小的尾翼斜置角时，使得上尾翼的两侧翼面再次处于不同的气流中，产生较大的压力差。同时发现，上下尾翼产生的压力差方向相反，形成了加速导弹转速的滚转力矩。对于左右尾翼而言，攻角的存在使得其两侧翼面产生较大的压力差，但由于左右尾翼压力差方向相同，因此对滚转力矩贡献较小。对比图10(a)和图10(b)，随着转速的增加，上下尾翼的两侧翼面压力差逐渐变小，由压力差产生的滚转力矩变小。从下文的定量分析中也可以得到该结果。

通过数值模拟获得了自旋尾翼弹在各种飞行条件下的滚转特性数据，如图11所示。由图11(a)可知，滚转力矩系数的绝对值随着马赫数的增大呈现先增大后减小的趋势，并且在Ma=1.1时其滚转力矩系数的绝对值最大。并且滚转力矩系数的绝对值随着转速的增大显著增大；同时发现图11(a)中的滚转力矩系数均小于零，这是由于尾翼斜置角γ=0°，导弹无法获得使其主动滚转的力矩(即导转力矩为0)而只存在阻碍其旋转的滚转阻尼力矩，故而滚转力矩系数均为负值。

对比图11(a)和图11(b)发现，两者随马赫数的变化趋势一致：其绝对值均随马赫数的增大呈现先增大后减小趋势。两者的区别在于：图11(b)中的滚转力矩系数均为正值，这是由于尾翼斜置角γ=1°时，该导弹存在一个使其主动滚转的力矩(即导转力矩)并且导转力矩的值大于滚转阻尼力矩的值，所以其滚转力矩系数表现为正值，在Ma=1.4附近滚转力矩系数取得最大值；同时还发现在图11(b)中，导弹转速越大，其滚转力矩系数越小，这是由于转速越大，导弹的滚转阻尼力矩也就越大。

对比图11(b)、图11(c) 和图11(d)可知：在同一马赫速和转速条件下，尾翼斜置角越大，则滚转力矩系数越大。这是由于尾翼斜置角越大，尾翼两侧的压力差也越大，导弹产生的导转力矩越大。

通过数值模拟获得了自旋尾翼弹在各种飞行条件下的滚转特性参数，如图12给出了不同Ma下导弹滚转力矩系数Cl随无量纲转速Ω·的变化趋势。

从图12(a)可知，当尾翼斜置角γ=0°时，导弹滚转力矩系数的绝对值随转速的增大呈现增大趋势，并与无量纲转速存在线性关系。对于同一转速条件下，在Ma=1.1下导弹的滚转力矩系数(即滚转阻尼力矩)的绝对值最大，Ma=2.0下导弹的滚转力矩系数的绝对值最小。同时Ma=1.1和Ma=0.9条件下的滚转力矩系数曲线基本平行。对比图12(b)和图12(c)、图12(d)可知：其滚转力矩系数均随无量纲转速的增大呈现减小趋势。

通过Matlab软件的拟合函数[b,bint,r,rint,stats]=regress(x,X)对导弹滚转数据进行拟合，得到了不同马赫数下，滚转力矩系数Cl与无量纲转速Ω·以及尾翼斜置角γ的关系式。如表5所示，发现滚转力矩系数与无量纲转速Ω·以及尾翼斜置角γ呈线性关系。将拟合值与仿真值对比发现，两者的最大平均误差小于9%，因此认为该关系式较为准确地描述滚转力矩系数与转速以及斜置角三者之间的关系。

根据表5的拟合函数得到了导弹飞行的平衡转速(即滚转力矩系数为零)，同时与工程算法得到的平衡转速进行对比，结果如表6～表8所示。当尾翼斜置角γ=0°时，旋转导弹在飞行过程中只存在阻碍其旋转的滚转阻尼力矩，即该导弹处于减旋过程中直到转速为0，因此其平衡转速为0。如表6所示，当尾翼斜置角γ=1°时，不同马赫数下的导弹的平衡转速不同，并且马赫数越大，其平衡转速越大。对比表6和表7、表8可知：对于同一马赫数而言，尾翼斜置角越大，其平衡转速越大。将数值计算方法得到的平衡转速和工程算法[14-16]得到的平衡转速进行对比发现，两者的最大误差不超过19%，考虑到所有的工程算法本身具有一定的误差，因此认为本文中通过拟合得到的经验公式(即表5中的拟合表达式)真实有效并且可以对平衡转速进行估算。

表5 拟合表达式

表6 γ=1°时导弹平衡转速

表7 γ=2°时导弹平衡转速

表8 γ=4°时导弹平衡转速

4 结论

1) 随马赫数增加，平均滚转力矩系数绝对值呈现先增加后减小的变化趋势，与常规滚转力矩变化规律一致。

2) 当尾翼斜置角为0°时，尾翼不产生导转力矩，其滚转力矩系数均为负值即滚转阻尼力矩；随转速增加，滚转阻尼力矩增强，且滚转阻尼力矩呈现典型的线性变化规律；当尾翼斜置角为1°、2°、4°时，其滚转力矩均为正值，且随转速的增大而减小。

3) 用Matlab软件拟合出滚转力矩系数与无量纲转速以及尾翼斜置角之间的函数关系式，并且将拟合值与仿真值以及工程算法得到的计算值进行比较，三者吻合较好，表明本文得到的经验公式真实有效。由于该公式可以快速地求出相应计算条件下的滚转力矩系数，以及相应的平衡转速，因此具有较大的实用价值。