鸭舵转速对双旋弹追随稳定性影响研究

摘 要：发展了适合双旋弹虚拟飞行仿真的高保真计算流体力学和刚体动力学耦合平台， 以此为基础研究了旋转鸭舵对双旋弹弹道追随稳定性的影响。 为了准确刻画双旋弹前后体的差动旋转效应， 将滑移网格算法引入自研非结构混合网格流场数值模拟程序HUND3D。 通过对双旋弹进行定轴转动非定常模拟， 考察了不同鸭舵转速下双旋弹的流动特征与气动特性。 通过耦合求解非定常雷诺平均NS方程和七自由度刚体动力学方程， 实现了双旋弹不同弯度弹道的虚拟飞行仿真， 分析了弹道追随过程的动力学机理， 并结合气动特性分析结果研究了控制鸭舵转速改善弹道追随稳定性的策略。 研究结果表明： 鸭舵旋转所产生的气动干扰， 能够显著影响弹体的侧向力与偏航力矩。 通过控制前体鸭舵转速以产生有利于弹道追随的偏航力矩， 能够在一定程度上改善弹道追随稳定性。

关键词：  双旋弹； 虚拟飞行仿真； 计算流体力学； 滑移网格； 侧向力； 偏航力矩； 追随稳定性

中图分类号：  TJ765

文献标识码： A

文章编号：  1673-5048（2024）02-0071-08

DOI： 10.12132/ISSN.1673-5048.2024.0069

0 引  言

随着现代战争对高精度、 低成本制导弹药的需求不断提高， 通过在制式弹药上增加简易制导装置以提升命中精度的弹道修正弹应运而生。 鸭舵式双旋弹作为一种弹道修正弹， 不仅可以利用后体高速旋转产生的陀螺稳定性保持弹轴方向的相对稳定， 防止弹体翻转， 还能通过控制前体鸭舵旋转来提高弹轴跟随速度方向变化的能力， 改善弹丸追随稳定性［1］。 对高速旋转弹丸而言， 提升陀螺稳定性与改善追随稳定性之间通常存在相互博弈的辩证关系： 弹体转速过低可能导致陀螺稳定性不足， 无法保证飞行姿态稳定； 转速过高， 弹丸陀螺稳定性过强， 则可能导致追随稳定性不足， 使得弹丸飞行迎角过大， 影响射程， 甚至出现弹底着地等现象。 因此， 要实现双旋弹高精度、 智能化制导控制的作战效能， 需要在陀螺稳定性和追随稳定性之间寻求平衡。

国内外学者围绕双旋弹的弹道追随稳定性开展了诸多研究。 Murphy［2］采用一阶线性化方法推导了无控弹箭攻角运动方程， 并对旋转弹飞行稳定需要满足的稳定性条件开展了研究。 考虑到前体修正引信对飞行稳定性的影响， 揭涛等［3］基于六自由度外弹道模型建立了采用减旋装置的双旋弹七自由度外弹道模型， 研究了不同射角下双旋弹的弹道追随稳定性， 获得了满足追随稳定性要求的最大射角。 冯单翔［4］基于攻角运动方程理论， 推导了双旋弹动力平衡角公式。 通过与普通旋转稳定弹经验公式对比， 发现双旋弹在弹道顶点处动力平衡角大于普通旋转稳定弹， 双旋弹的追随稳定性较差。

Moss等人［5］曾在研究中指出， 当旋转弹自转速率接近偏航速率时， 会出现灾难性偏航等非线性现象。 对于鸭式布局双旋弹， 常江思等［6］基于小攻角条件推导出的攻角运动方程， 发现前体鸭舵转速与攻角圆运动频率相近时， 会出现共振等不利于飞行稳定的现象。 随着研究的深入， Yin等［7］通过对比URANS（Unsteady Reynolds- Averaged Navier-Stokes）方法与外弹道学理论公式下旋转弹的气动力系数， 发现在较高自旋速率下两种方法获得的气动力系数存在显著差异， 表明理论公式方法在预测气动力方面存在一定的局限性。 基于滑移网格技术与S-A（Spalart-Allmaras）湍流模型， Liu等人［8］研究了后体自转速度、 马赫数以及攻角对双旋弹气动特性的影响， 发现前体翼尖涡能够使得后体气动特性发生周期性变化， 且变化频率等于前体自旋速率。 黄智康等人［9］采用滑移网格技术对超声速状态下双旋火箭弹进行了数值模拟， 分析了鸭舵对弹体侧向力的影响。 研究发现， 鸭舵组件显著影响弹体侧向力， 且加装鸭舵后弹体侧向力系数的非线性随攻角的增加而变强。

综上所述， 目前有关弹道追随稳定性研究大多采用小迎角线性化和准定常假设等简化方法进行研究， 忽略了飞行过程中气动力和力矩的非定常与非线性。 近年来， 数值虚拟飞行仿真［10-11］（Numerical Virtual Flight， NVF）已在亚声速旋转弹［12］、 超声速旋转弹［13-14］等高速自转飞行体气动力及飞行历程预测中得到广泛应用， 并取得满意结果。 为研究复杂气动力下双旋弹的追随稳定性， 有必要通过非定常CFD驱动的数值虚拟飞行仿真方法， 精确预测双旋弹所受非定常气动力， 高保真地获得飞行过程中的动力学响应， 进而更准确地进行双旋弹飞行稳定性评估， 研究旋转鸭舵气动干扰对弹道追随稳定性的影响。 为此， 本文建立了适用于双旋弹飞行稳定性研究的高保真数值虚拟飞行仿真平台。 通过对双旋弹进行非定常定轴转动数值模拟以及不同弯度弹道下虚拟飞行仿真， 分析不同前体转速下的气动特性与追随过程动力学机理。 在此基础上， 提出采用鸭舵反转来改善弹道追随稳定性的控制策略， 并通过仿真验证进行效果评估。

1 计算方法

本文通过耦合求解非定常流动控制方程（CFD模块）与七自由度刚体动力学方程（RBD模块）， 实现双旋弹数值虚拟飞行仿真。

1.1 非定常流动控制方程

采用自主研发的MPI并行非结构混合网格N-S方程求解程序HUNS3D［15］ ， 对双旋弹飞行过程中的非定常绕流进行数值模拟， 进而获得非定常气动力。 为了考虑多自由度运动控制体网格相对惯性坐标系的运动， HUNS3D采用任意拉格朗日-欧拉方法（Arbitrary Lagrange-Euler， ALE）描述的积分形式流动控制方程为

tΩQdV+ΩF（Q，Vgrid）·ndS=ΩG（Q）·ndS（1）

式中： Ω代表控制体； Ω代表控制体边界； 守恒变量Q=ρρuρvρwρET， ρ为流体密度， u、 v、 w分别为流体x、 y、 z方向的速度， E为控制体的总内能； Vgrid为网格运动速度； n为控制体表面外法线单位向量； 无粘通量项F（Q，Vgrid）·ndS与粘性通量G（Q）·ndS的具体形式见文献［15］。 网格运动速度Vgrid采用二阶向后差分公式［16］计算获得。

对流动控制方程式（1）中的无粘通量F（Q，Vgrid）·ndS与粘性通量G（Q）·ndS， 使用格心型有限体积法进行离散逼近［17］， 可以得到有限体积单元上半离散形式N-S方程：

Qit=－1Ωi（Ri－Rvi）（2）

其中， Ri与Rvi分别为当前单元i的无粘通量项与粘性通量项对残值的贡献。 通过对空间离散后的（2）式采用全隐式双时间格式［18］以完成时间推进， 具体可参阅文献［19］。

1.2 滑移网格技术

通常， 双旋弹由于前后体转速不同存在物面相对运动的情况。 为准确模拟前体鸭舵相对后体旋转对流动的影响效应， 本文在CFD求解程序HUNS3D中发展了滑移网格技术。 考虑到CFD计算过程中需要完成各面单元的通量计算， 而滑移边界面单元只具有其所属网格内侧的流场值， 必须与其他进程进行通信以获得单元外侧流场值。 为实现滑移边界两侧计算数据的通信， 本文采用图1所示的滑移网格算法策略。 以面单元La2为例， 将滑移边界另一侧距离最近的面单元Rb2作为贡献单元， 以贡献单元Rb2及其所有相邻面单元流场值作为插值基， 插值获得La2滑移面外侧的流场值， 从而进行后续通量计算。

为验证所开发的滑移网格算法， 本文进行了近地悬停共轴四旋翼数值模拟。 基本外形参数如下： 旋翼截面翼型为NPL9615， 弦长c=0.18 m， 展长0.625 m， 根梢比为1， 旋翼半径R=1.105 m， 桨距为15°， 无扭矩， 四片旋翼旋转对称。 计算状态如下： 旋翼离地高度h/R=0.52， 翼尖马赫数Ma为0.56， 雷诺数Re为1.44×107， 参考长度取旋翼半径R。 计算网格如图2所示， 图中白色与红色网格分别为滑移边界与旋翼表面网格， 旋翼附面层第一层网格高度Y+为0.8， 法向高度增长率为1.15， 旋翼尾迹区的网格尺度为0.09c。

近地悬停共轴四旋翼滑移边界处流场如图3所示。 图3（b）旋翼翼尖涡在滑移网格边界两侧光滑连续， 以及图3（c）径向截面中流场马赫数等值线的连续分布， 表明了本文所使用的滑移网格数据传递技术能够正确实现滑移边界两侧流场值通信。 图4展示了近地悬停共轴四旋翼数值计算与试验的翼尖涡轨迹空间位置对比图。 图中， 数值计算所得翼尖涡轨迹与试验测量结果［20］基本一致。 此外， 数值计算所得旋翼升力系数0.092 3与试验测量值0.090的差异很小， 表明本文计算结果是准确的。

1.3 七自由度刚体动力学方程

为建立七自由度刚体动力学方程， 本文引入了有关坐标系以描述双旋弹运动与动力学参数。 图5为双旋弹地面坐标系OgXgYgZg、 前体体轴系ObXbYbZb以及弹道坐标系OVXVYVZV的示意图。 地面坐标系原点Og位于双旋弹飞行起始点， Xg轴沿水平线指向初始速度方向， Zg轴铅垂向下， Yg轴由右手定则确定。 前体体轴系原点Ob固联于双旋弹全弹质心， Xb轴沿弹体轴线方向指向前方， Yb轴与Zb轴取为鸭舵面平面的法向量方向。 弹道坐标系原点OV与前体体轴系原点Ob重合， XV轴沿质心平动速度方向指向前方； ZV轴位于铅垂平面内， 垂直XV轴且方向向下； YV轴通过右手法则确定。 通过前体体轴系ObXbYbZb与双旋弹地面坐标系OgXgYgZg的相对方位， 可以给出双旋弹弹轴相对地面的飞行姿态角定义： 俯仰角θ为Xb轴与水平面OgXgYg的夹角， 偏航角ψ为Xb轴在水平面OgXgYg的投影与Xg轴的夹角。 类似地， 可以获得质心速度相对地面的方位θV与ψV的定义。 基于上述坐标系与姿态定义， 可以计算得到外弹道学理论［1］中复攻角： δ=α+iβ=（θ－θV）+i（ψ－ψV）。

本文基于前体体轴系ObXbYbZb将双旋弹前体与后体两组六自由度动力学方程简化为一组七自由度动力学方程。 在前体体轴系下， 质心平动动力学方程矢量形式为

Fb=mdVbdt+ω×Vb（3）

式中： m为双旋弹全弹质量； Fb为质心在前体体轴系下合外力矢量； Vb为质心在前体体轴系下平动速度； ω为前体体轴系相对地面坐标系的旋转角速度。

由于双旋弹前后体间存在约束， 因此将旋转自由度最终简化为俯仰自由度、 偏航自由度以及前后体的两个滚转自由度。 双旋弹绕全弹质心转动的动力学方程具体形式如下：

Mfx=Ifxdωfxdt

Max=Iaxdωaxdt

My=Iydωydt+Ixω-xωz－Izωfxωz

Mz=Izdωzdt+Iyωfxωy－Ixω-xωy

ω-x=Ifxωfx+IaxωaxIx

（4）

式中： M为相对全弹质心的合力矩； I为转动惯量； 下标f与a分别代表双旋弹前体与后体。 对于前体受电机控制维持在恒定转速的情况， 本文将式中前体所受力矩Mfx叠加至后体Max后， 对Mfx进行置零操作。

为求解式（3）与式（4）所组成的常微分方程组， 本文采用了一种改进的Adams预估校正法［16］： 通过已知时间步气动力插值获得校正步所需的下一时间步气动力， 以提高计算效率。

1.4 CFD/RBD耦合流程

图6为数值虚拟飞行仿真平台中CFD/RBD耦合计算流程图。 CFD模块在完成当前时间步计算后输出气动力与力矩系数； RBD模块将CFD输出有量纲化后作为输入参数， 求解获得下一时间步的位移与姿态角变化。 根据RBD求解所得动力学参数， CFD模块对计算网格做出相应的刚体变换， 并计算网格运动速度。 CFD和RBD按照上述流程进行耦合迭代， 直至满足结束条件。