伴飞弹是由载体弹射并具有一定伴飞能力的飞行器[1]，为了能够在短时间内与载体分离并保持一定的相对距离飞行，伴飞弹通常垂直于载体飞行方向弹射脱离载体并通过自身气动力[2]、推进发动机[3]或被动拖曳力[4]来迅速调整到稳定飞行姿态，其中采用自身气动力调整的工作方式要求弹体具有较高的自稳定性能，而设计满足战术要求的弹体结构就需要对弹体在大攻角变化范围内的气动特性进行研究。
飞行器大攻角飞行状态下的气动力计算与绕流场仿真问题非常复杂，利用传统的工程计算方法已经不能得到准确的结果，目前国内外研究人员主要通过风洞实验和数值仿真的方法对这一问题进行研究，并得到了一定的研究成果。Fairlie[5]最早研究了头部为半球形的细长旋成体在大攻角情况下的绕流情况，发现当攻角大于某值时旋成体头部会产生一对对称的角涡；Tobak和Peake[6]基于奇点理论[7]对大攻角旋成体绕流场的三维拓扑结构进行了研究，得到了旋成体背风面气流分离再附着的典型结构；Clainche[8]则研究了不同雷诺数和攻角下半球-圆柱体背风面绕流场的形成机理；随后人们发现攻角继续增大后，细长旋成体的绕流逐渐发展为非对称涡结构[9]，且非对称的绕流结构对旋成体产生了不断变化的侧向力[10]；为了研究侧向力的变化规律，Zhu[11]和Qi[12]分别对尖头-圆柱体和钝头-圆柱体进行了试验和数值仿真研究，得到了头部微小突出物对侧向力的影响规律。从现有的研究成果来看，人们对细长旋成体这样结构较为简单的飞行器在大攻角状态下的绕流场和受力问题已经有了较为详细的研究，而对于带有尾翼的结构更为复杂的飞行器，其在大攻角状态下的绕流场和受力情况则会产生很多变化，仍需要进行详细研究。
本文研究对象为某外形复杂、长径比较小的无控伴飞弹，由于其尾翼结构的非对称性，导致其气动特性与常规尾翼弹有所差别，用传统的工程计算方法不能精确地计算出其所有气动特性，因此采用数值仿真的方法对其在0°～90°攻角范围内进行计算，得到了各个攻角下的气动参数和绕流场分布情况，并发现该弹体结构在一定攻角范围内的静稳定性能较差，随后提出了进一步的改进方案，为工程设计和实际应用提供了理论和数据支持。
1 伴飞弹全六面体结构网络
本文研究的伴飞弹头部为截卵形，弹身为圆柱形，弹径为D=47.2 mm，弹体尾部有四片折叠尾翼及相应的折叠机构，其中四片尾翼关于弹底面中心点对称，在弹头与弹身结合处设置了两道密封槽用以密封发射气体，具体尺寸如图1所示。
图1 伴飞弹示意图
由于结构网格相较于非结构网格具有高精度、高效率、贴体性好、各向异性好等优点，更适用于计算本文所研究的大攻角飞行弹丸复杂绕流场，因此本文采用全六面体结构网格。由于要计算伴飞弹在大攻角(αmax=90°)飞行时的流场情况，为了避免边界对计算结果的影响，所以将外流场范围设置的较大，为伴飞弹直径的30倍，如图2所示；采用O形结构网格包围弹身和尾翼，对压力变化剧烈的部位(如弹头部、密封槽、尾翼)进行了加密；为确保Y+≤1，壁面法向第一层网格高度为2×10-5，边界层增长率为1.2；通过精细调整网格拓扑结构和网格分部，得到了质量较好的网格，网格总数目为8 940 000。
图2 计算网格分布示意图
2 数值方法
利用有限体积法求解三维定常雷诺平均N-S方程组(RANS)，在笛卡尔坐标系中，其形式为：
(1)
式中：Q为守恒变量； F、G和Η为对流项；Fv、Gv和Hv为黏性项； t为时间；x、y和z为坐标系的3个方向。
式(1)为控制方程，采用二阶迎风格式对动量与时间进行离散，基于SIMPLE算法进行压力速度耦合。
湍流模型采用SST k-ω模型，对于近壁面区域和远离壁面的湍流区域均有较好的适应能力，能够准确地计算逆压梯度的流体分离情况，具体方程如下：
(2)
(3)
方程中具体参数的含义见文献[13]。
进出口边界条件选用压力远场条件，来流马赫数为0.6，压力为101 325 Pa，温度为300 K，流体选用理想气体。壁面边界采用等温壁面假设和无滑移条件。
对文献[14]所述长径比为4的试验模型进行了数值计算，计算结果与试验所测结果如图3所示。可以看出计算所得俯仰力矩系数与试验数据具有较好的一致性，平均误差小于4.3%；计算所得侧向力系数相对试验数据存在一定的误差，这是由于试验模型存在微小缺陷或加工误差，与理想模型并不完全一致，而这些微小误差会对绕流场产生较大的影响，导致计算所得侧向力系数与试验数据相差较大，但侧向力系数曲线整体变化趋势和试验数据基本一致。
图3 数值计算结果与文献[14]试验结果
3 仿真结果
通过上述数值方法计算了伴飞弹在0～90°攻角范围内(每5°为一个计算状态)的绕流场及气动力参数，全面地展示了伴飞弹在各个攻角下的受力情况：
3.1 中垂面速度与压力云图
图4为攻角为0°～90°的弹体表面压力及绕流场云图，其中对称面为马赫数云图，弹身表面为压力云图。
图4(a)是攻角为0°时的弹体表面压力及绕流场云图，可以看出弹体绕流场关于弹轴对称。弹头前部由于流动滞止的原因产生了局部高压，气体绕过弹头截平面后在卵形弹头部分逐渐加速，相对应的弹头表面压力也逐渐降低。直到接近第一道密封环时，由于绕流受阻，弹体表面压力有所上升，密封环后的压力则逐渐降低直至到达弹身圆柱段，而弹身圆柱段表面压力则保持均匀分布。到达尾翼前缘时气体被压缩，弹体表面压力升高，而越过尾翼前缘后的气体迅速膨胀增速，随后在翼尖涡、翼根处复杂涡流的作用下，尾翼表面压力呈现出不均匀的分布。最后在弹体尾部产生了关于轴线对称的尾涡。
图4(b)是攻角为15°时的弹体表面压力及绕流场云图，可以看出流场结构不再关于弹轴对称，弹头部的高压区向迎风面偏移，导致弹头背风面产生低压区，而且弹头背风面气体流速相较于零攻角状态也有所增大，但未出现明显的流动分离现象。在密封环的迎风面仍然产生了和零攻角时相似的流动状态，而密封环背风面局部气体流速则未出现明显的变化。弹身圆柱段则呈现出两侧面压力低于迎风面与背风面压力的情况。在尾翼段，由于迎风面气体受到压缩，导致表面压力高于背风面，另外尾翼折叠槽表面的压力呈现出不均匀的变化，这是由于气体在槽内产生涡流所引起，具体分析将在下文中指出。弹体尾部的流场则由于在小攻角情况下弹底面靠近背风面侧真空度大于靠近迎风面侧的原因，导致背风面气流在底部膨胀偏转角度更大，呈现出“下垂”而非“上翘”的不对称状态。
图4(c)是攻角为30°时的弹体表面压力及绕流场云图，可以看出弹头部高压区进一步向迎风面偏移，背风面低压区范围扩大，且绕流产生了明显的分离现象。弹身圆柱段两侧表面压力相较于15°攻角时更低，而迎风面表面压力则高于背风面。在尾翼段迎风面与背风面表面压力差别更明显，且上垂尾两侧表面压力相较于15°攻角时也降低很多，基本与两侧水平尾翼的背风面压力相同。弹体尾部流场相较于15°攻角时也产生了较大变化，流场结构开始呈现出“上翘”趋势，这说明随着攻角的增大，弹底面靠近迎风面侧气流的径向速度分量开始大于背风面的值，导致流场呈现图中所示结构。
图4(d)是攻角为45°时的弹体表面压力及绕流场云图，可以看出弹头部高压滞止区已全部移动到弹头迎风面，气流分离区也进一步扩大，在轴向距离弹头截面约28.5 mm的位置气流再次附着。随着攻角的不断增大，弹体迎风面压力不断升高，而弹体两侧面和背风面压力均进一步降低。在尾翼段，由于气体涡流的作用，上垂尾左侧面(弹头方向为前方)前部的压力相较于其他情况要低很多，这将产生额外的侧向力和滚转力矩，影响到弹体的姿态，具体分析见下文。
图4(e)～(g)分别是攻角为60°、75°和90°时的弹体表面压力及绕流场云图，可以看出随着攻角的不断增大，弹体表面压力按照上文所述的趋势在不断变化，弹体下游的低速流场区域也由于背涡逐渐远离弹体而不断变大。值得指出的是，在攻角为75°时，弹头前部的气体流速最大，导致弹头背风面低压区范围也最大。
图4 弹体表面压力及绕流场云图
3.2 背风面涡流图
由于随着攻角的增大，弹体背风面的绕流会由对称结构转变为非对称结构，使弹体受到侧向力，对弹体的稳定性产生一定影响，因此需要对背风面的绕流情况进行分析。
图5(a)～(d)为弹体周围涡量云图，展示了在各个攻角时弹体不同位置的绕流结构，图中所示位置1在x/L=0.23处(L为弹体总长)；位置2在x/L=0.451处；位置3在x/L=0.674处；位置4在x/L=0.897处。
图5 不同攻角时弹体周围涡量云图
由图5(a)可以看出，当α=30°时，弹身部分(位置1、2、3)背风面产生了脱离壁面的对称背涡，由于尾翼结构的不对称导致弹翼部分(位置4)的背涡是非对称的，因此会在尾翼段产生部分侧向力和力矩；由图5(b)可以看出，随着攻角的增大，弹体背风面的对称背涡逐渐远离弹体(其中位置3处的背涡已经完全脱离了弹体)，尾翼段的非对称背涡也产生同样的变化；由图5(c)可以看出，当攻角增大到45°时，弹身背风面原来基本对称的背涡也发展成了非对称的，位置1处左侧(由弹尾朝弹头方向看)的背涡有一部分被右侧的涡压住，而在位置2和3处涡的非对称性逐渐增大，弹体所受侧向力和力矩也相应增大；随着攻角的继续增大，弹体背风面的背涡逐渐脱离弹体，由图5(d)可以看出，当攻角增大到75°时，弹体背风面(位置2、3、4)的背涡已经完全从壁面脱离，且非对称性极强。
3.3 气动力与力矩：
通过流场速度、涡量和弹体表面压力云图可知伴飞弹在大攻角范围内所受气动力呈非线性变化，本节通过数值计算得到伴飞弹在Ma=0.6下主要气动参数并进行了详细分析。计算结果以气动力系数和气动力矩系数的形式表示，其中升力系数与阻力系数是基于速度坐标系(来流方向为x轴正向，垂直于来流方向向上为y轴正向，z轴正向由右手定则确定，原点位于弹体质心处)，其余参数是基于弹体坐标系(沿弹轴指向弹头方向为x轴正向，垂直于弹轴向上为y轴正向，z轴正向由右手定则确定)，如图6所示。
图6 弹体主要气动参数随攻角变化曲线
由图6(a)中升力系数曲线可以看出伴飞弹所受升力随攻角先增大后减小，在攻角为50°时达到最大值。而由阻力系数曲线可以看出弹体所受阻力则随攻角的增大不断增大，直到攻角为75°左右时达到最大值，随后阻力有所减小，这是由于在攻角为75°左右时弹体背风面的低压区范围最大，使得迎风面与背风面所受合力相差最大，导致阻力达到峰值，当攻角继续增大后，弹体背风面绕流速度相对降低，使得表面压力有所增大，因此迎风面与背风面所受合力差减小，阻力也相应减小。
图6(b)所示为弹体侧向力系数曲线，可以看出在攻角逐渐增大到75°之前，侧向力的大小和方向随着攻角的增大呈现出周期性的变化，大约以25°为一个周期。在攻角小于30°时，弹身部分绕流关于弹体中垂面对称分布，而由于弹翼部分不是关于中垂面对称的，随着攻角的增大，这部分的绕流则不断变化，导致侧向力的指向呈周期性的变化。在攻角大于30°后，弹身部分背风面的绕流变为非对称的，弹身表面左右两侧的压力差导致侧向力变化幅度增大，但是方向仍延续之前的周期性变化，说明弹体周围的绕流场结构总体呈现出周期性的变化。直到攻角大于75°之后，侧向力的大小不断波动，方向则始终指向左侧(z轴负向)，这是由于在攻角大于75°后，弹身背风面的涡流不再是定常的，因此所产生的侧向力不断波动，而由于弹翼部分迎风面结构的不对称性，导致侧向力总体偏向左侧。
图6(c)所示为弹体所受俯仰力矩系数变化曲线，可以看出在攻角从0°～90°的变化范围内，只有在攻角很小(α<5°)时弹体会受到正的“仰头”力矩，当攻角大于5°后弹体始终受到负的“低头”力矩，说明弹体是静稳定的。在攻角小于20°时，弹体所受俯仰力矩随攻角的增大而减小；而当攻角在20°～50°范围内，弹体所受俯仰力矩的大小随攻角的增大而发生振荡，但方向仍保持不变；当攻角大于50°时，俯仰力矩的值又单调减小。俯仰力矩系数之所以产生上述变化规律是由于当攻角小于20°时弹体背风面未产生背涡，此时水平尾翼背风面的表面压力随着攻角的增大而逐渐变小，迎风面的表面压力则逐渐增大；当攻角在20°～50°范围内时，弹体背风面产生了非稳态的背涡，这些背涡作用于水平尾翼背风面，导致其表面压力产生了非线性的增大，与迎风面的压力相抵；攻角继续增大后，弹体背风面产生的背涡逐渐脱离弹体，又使得尾翼背风面表面压力迅速减小，弹体背风面流场的上述变化最终造成了图6(c)所示俯仰力矩系数的变化规律。另外，虽然弹体在整个攻角变化范围内是静稳定的，但是在相当大的攻角范围内(25°<α<60°)弹体的俯仰力矩系数的绝对值都小于0.2，甚至在55°攻角时接近于零，这不利于伴飞弹姿态迅速恢复稳定。因此有必要对弹体结构进行改进，采用质心位置前移或增大尾翼面积以便将压心位置后移等方式，提高伴飞弹的静稳定性。
图6(d)所示为弹体滚转力矩系数曲线，可以看出在攻角小于5°时滚转力矩系数大于0，说明弹体受到顺时针(由弹尾朝弹头方向看)的力矩，在其余攻角变化范围内滚转力矩系数均小于0，说明弹体始终受到逆时针方向的滚转力矩。滚转力矩主要由弹体所受升力诱导产生[15]，由于左右两侧的水平尾翼的不对称性，导致两翼所受升力不均，从而产生滚转力矩。而上垂尾左右两侧面的压力差也会对滚转力矩产生一定影响，如上文所述，当攻角为45°左右时上垂尾左侧面出现低压区，导致逆时针的滚转力矩明显增大，图6中滚转力矩系数曲线也反映出这一点。
3.4 尾翼结构改进：
为了提高伴飞弹的静稳定性，综合考虑弹体结构后，提出改进两种方案：1) 将尾翼弦长增大5 mm(16.7%)；2) 将翼展增大10 mm(11.8%)。针对这两种方案进行仿真计算后，得到了各个攻角下的气动力参数，如图7所示。
图7 尾翼改进后的气动力系数曲线
图7(a)为各个弹型的升力、阻力系数随攻角变化曲线对比图，可以看出两种改进方案均使得伴飞弹的升力、阻力系数有所增大。其中，升力的增大量首先随着攻角的增大而增大，当攻角大于45°后，增大量则随着攻角的增大而减小；阻力系数的增大量则随着弹体攻角的增大而持续增大。
图7(b)展示了各个弹型的俯仰力矩系数随攻角变化的曲线对比，可以看出尾翼改进后，伴飞弹俯仰力矩系数的绝对值在各个攻角下均有所增大，而且方案二的增量更大，说明增大尾翼的翼展能够更加有效的提高伴飞弹的静稳定性。
为了能够进一步对两种改进方案进行定量分析，引入一个参数：系数总增量ΔC，具体表达式如下：
(1)
式中: C(α)为改进方案在攻角为α时的气动系数； C0(α)为原型弹在攻角为α时的气动系数。表1列出了两种改进方案升、阻力系数和俯仰力矩系数的系数总增量，从表中可以看出两种改进方案的俯仰力矩系数的系数总增量均明显大于升、阻力系数的系数总增量，说明这两种改进方案对弹体静稳定性的提升效果要明显高于对升、阻力系数的影响；另外，从表中还可以看出翼展增大方案对弹体静稳定性的提升也有优于弦长增大方案，与前文所述结论一致。
表1 两种改进方案主要气动参数的系数总增量
CxCyCmz弦长增大0.633.8864.95翼展增大0.882.5580.37
4 结论
1) 弹体头部绕流在攻角为30°时表现明显的流动分离，分离区随着攻角的增大不断扩大甚至完全扩大到整个弹体背风面，导致流动脱离壁面。在攻角为75°时弹头背风面低压区范围最大，弹体所受阻力也达到最大；
2) 弹体背风面绕流在攻角为35°时为非对称，随着攻角的增大背风面非对称绕流呈现周期性变化，导致弹体所受侧向力也出现以25°为周期的变化，直至攻角大于75°；
3) 在攻角大于5°后，弹体始终保持方向不变，说明弹体是静稳定的，但在很大攻角变化范围内俯仰力矩系数的绝对值较小，不利于伴飞弹迅速调整到稳定状态，应调整质心位置、改变气动力压心位置等方法保持稳定状态。
4) 增大尾翼翼展比增大尾翼弦长更能有效增加伴飞弹的静稳定性。
[1] 王庆华，田桂荣，李军.红外伴飞诱饵研究初探[C]//第十二届光电对抗与无源干扰学术交流会论文集.北京：中国电子学会，2002：99-102.
[2] HERBAGE D W，HEINSOHN W H.Decoy flare[P].U.S.Patents：5400712.1995.
[3] TOWNING J N，SUTLIEFF T D，PELHAM P G.Propelled pyrotechnic decoy flare[P].U.S.Patents：5561260.1996.
[4] 王鹏.拖曳式诱饵弹释放过程绳索动力学仿真研究[D].长沙：国防科学技术大学，2011.
[5] FAIRLIE B.Flow separation on bodies of revolution at incidence[C]//Proceedings of the 7th Australian Hydraulics and Fluid Mechanics Conference.Brisbane，Australia：Institution of Engineers，1980：18-22.
[6] TOBAK M，PEAKE D.Topology of two-dimensional and three-dimensional separated flows[C]//Proceedings of the 12th Fluid and Plasma Dynamic Conference.Virginia，USA：AIAA，1979：1-29.
[7] DALLMANN U.Topological structures of three-dimensional vortex flow separation[C]//Proceedings of the 16th Fluid and Plasma Dynamics Conference.Danvers，USA：AIAA，1983：1-10.
[8] CLAINCHE S L，RODR
GUEZ D，THEOFILIS V，et al.Formation of three-dimensional structures in the hemisphere-cylinder[J].Aiaa Journal，2016，54(12)：1-11.
[9] DEGANI D.Numerical investigation of the origin of vortex asymmetry[C]//Proceedings of the 28th Aerospace Sciences Meeting.Moffett Field，USA：AIAA，1990：1-12.
[10] 曾广存，秦伊贤，丁庆国，等.大攻角非对称涡特性及消除侧向力和力矩方法[J].空气动力学学报，1988，6(1)：137-141.
[11] ZHU Y，YUAN H，LEE C.Experimental investigations of the initial growth of flow asymmetries over a slender body of revolution at high angles of attack[J].Physics of Fluids.2015，27(8)：1-13.
[12] QI Z，WANG Y，BAI H，et al.Effects of micro perturbations on the asymmetric vortices over a blunt-nose slender body at a high angle of attack[J].European Journal of Mechanics / B Fluids.2018，68(2018)：211-218.
[13] MENTER F R.Two-equation eddy-viscosity turbulence models for engineering applications[J].Aiaa Journal,1994，32(8)：1598-1605.
[14] LOPERA J，NG T T，PATEL M P，et al.Forebody geometry effects on the flow of a blunt-nose projectile at high alpha[J].Journal of Aircraft,2007(44)：1906-1922.
[15] 郑健，周长省，张磊.弧形翼-身组合体零攻角流场的数值仿真研究[J].系统仿真学报，2008，20(13)：3401-3403，3419.