引言

飞行器在大气环境中真实飞行时，其运动是多个自由度的耦合.传统上研究多自由度耦合运动效应时，一般是把单自由度运动结果线性叠加起来^[1].但随着现代飞行器机动性和敏捷性的提升，特别是大攻角飞行时，这种处理方式是否可行仍存在疑问.如Wang等^[2]的研究发现，采用常规准定常模型和单自由度非线性叠加模型进行飞机横向机动运动仿真时结果与试飞结果有很大差别.贾区耀等^[3-4]在对飞行器大气自由飞的结果和风洞实验数据进行对比分析时发现，风洞实验或数值计算表明是稳定的某些飞行器，自由飞的结果却是不稳定的，分析认为自由度不相似可能是导致这些差异的重要原因之一^[5].

在飞行器的运动多自由度耦合效应研究中，纵横向耦合效应是研究的重要方向之一.在实验研究方面，国内如南京航空航天大学黄达等^[6-8]在低速风洞中研究了三角翼外形俯仰-滚转耦合特性，发现机动飞行时多自由度运动的飞行器气动特性比单自由度运动时要复杂的多，运动耦合时的气动特性与两个单自由度运动时的气动特性的叠加结果相比有一定差别.唐敏中等^[9]分析了80$^\circ$三角翼不同相位角和不同频率对耦合气动特性的影响，发现双自由度的耦合运动比单自由度运动具有更明显的动态迟滞特性, 会进一步推迟分离涡的破碎, 推迟开始滚摆的攻角.近年来，随着多自由度运动耦合效应逐步受到重视，国内多家单位发展了相应的风洞试验技术，具备了对相关问题的研究能力^[10-15].

在数值模拟方面，杨小亮等^[16-18]针对三角翼外形采用数值模拟手段研究了多种因素对俯仰/滚转运动特性的影响，发现俯仰耦合滚转运动会导致飞行器升力降低，横向稳定性变差.郭迪龙等^[19]数值模拟了76$^\circ$三角翼强迫俯仰/强迫滚转耦合运动过程，发现俯仰滚转耦合运动时，三角翼上表面的涡分布的非对称性将产生横侧方向的偏航力矩和滚转力矩，滚转力矩和偏航力矩随着滚转振幅角和滚转缩减频率的增大而增大，但对法向力影响不大.

国外也有多家机构和研究人员采用实验^[20-23]或计算^[24-27]的手段研究了飞行器的多自由度运动耦合效应.总体来看，由于俯仰、滚转耦合运动的复杂性，虽然取得了一定的研究成果，但仍需要开展更深入的研究工作.另外，研究的对象多集中在三角翼外形，对导弹类外形的研究结果较少，而运动耦合效应对气动外形的依赖度很大.本文通过建立刚体六自由度动力学方程与Navier-Stokes (N-S)方程的耦合求解技术，研究滚转静稳定和静不稳定时，方形截面弹在不同振动频率下的俯仰振荡过程对滚转静、动态稳定性的影响，在此基础上，开展俯仰振荡过程对滚转运动特性的影响研究.

1 数值方法和计算模型 1.1 数值方法

在数值求解非定常Navier-Stokes方程时，空间离散格式采用原始变量NND (none oscillation，none free parameter, dissipative difference scheme)格式^[28]，限制器选用min-mod限制器；非定常时间推进采用Jameson的双时间步方法^[29]；在求解刚体动力学方程时，采用显式四阶龙格-库塔方法进行求解.采用松耦合的方式解决N-S方程和六自由度刚体动力学方程的耦合求解问题，两组方程分别独立求解，在时间域上交错推进，从而获得耦合系统的响应.

1.2 计算模型与来流条件

计算模型为方形截面弹外形^[30]，由一个圆弧形头部和方形截面的弹身组成，参考图 1.弹身总长$L=13D$，$D=93.98$ mm，头部为切线圆弧，长3$D$，圆弧半径$R=6.717D$.弹身截面为1$D\times $1$D$的正方形.弹身尾部有四片"十"字布局的三角小翼，小翼的长和高均为$D$，宽为0.1$D$.

计算马赫数$Ma=2.5$，以弹体全长为参考长度的雷诺数$Re=1.6\times $10$^{7}$.力矩参考点在弹长的0.47处.飞行器的转动惯量取$J_{x}$=0.1，$J_{y}$=2.0，$J_{z}$=2.0.

1.3 刚体动力学方程分析

首先定义体轴系.体轴系固连于飞行器，原点$O$为飞行器的质心；$Ox_{1}$轴与飞行器纵轴重合，指向头部为正；$Oy_{1}$轴在纵向对称面内，垂直于$Ox_{1}$轴且向上为正；$Oz_{1}$轴按右手系确定.

由于飞行器一般是关于$z_{1}$=0平面对称的，而对此类面对称的飞行器，飞行器绕质心转动的动力学方程可写为如下标量形式

$\left.\begin{array}{lll} \dfrac{{\rm d}\omega _x }{{\rm d}t} = a_1 \omega _x \omega _z + a_2 \omega _y \omega _z + \bar {M}_{x1} \\ \dfrac{{\rm d}\omega _y }{{\rm d}t} = b_1 \omega _x \omega _z + b_2 \omega _y \omega _z + \bar {M}_{y1} \\ \dfrac{{\rm d}\omega _z }{{\rm d}t} = c_1 \omega _x \omega _y + c_2 \omega _x^2 + c_3 \omega _y^2 + \bar{M}_{z1} \\ \end{array}\right\}$

(1)

其中，

$\left.\begin{array}{lll} a_1 \equiv \dfrac{I_{xy} \left( {I_z-I_x-I_y } \right)}{I_x I_y-I_{xy}^2 }, ~&a_2 \equiv \dfrac{I_{xy}^2-I_y I_z + I_y^{2} }{I_x I_y-I_{xy}^2 }, \\ b_1 \equiv \dfrac{I_x I_z-I_x^2-I_{xy}^2}{I_x I_y-I_{xy}^2 }, &b_2 \equiv \dfrac{I_{xy} \left( {I_x + I_y-I_z } \right)}{I_x I_y- I_{xy}^2 }, \\ c_1 \equiv \dfrac{I_x-I_y }{I_z}, c_2 \equiv \dfrac{I_{xy}}{I_z }, &c_3 \equiv-\dfrac{I_{xy}}{I_z}, \\ \bar {M}_{x1} \equiv \dfrac{I_y M_{x1} + I_{xy} M_{y1}}{I_x I_y-I_{xy}^2 }, &\bar {M}_{y1} \equiv \dfrac{I_{xy} M_{x1} + I_x M_{y1}} {I_x I_y-I_{xy}^2}, \\ \bar{M}_{z1} \equiv \dfrac{M_{z1}}{I_z} \\ \end{array}\right.$

式中，$\omega_{x}$，$\omega _{y}$，$\omega _{z}$为转动角速度的分量，$I_{x}$，$I_{y}$，$I_{z}$，$I_{xy}$为惯性张量的分量，$M_{x1}$，$M_{y1}$，$M_{z1}$为力矩在体轴系下的分量.对选择的方形截面导弹外形或大多数圆形截面导弹外形而言，外形在关于$z_{1}$=0平面对称的同时，也关于$y_{1}$=0平面对称；因此除$I_{xz}$和$I_{yz}$为零外，$I_{xy}$亦为零.再结合$I_{y}=I_{z}$，则式(1)可化简为

$\left.\begin{array}{l} \dfrac{{\rm d}\omega _x }{{\rm d}t} = \dfrac{M_{x1} }{I_x} \\ \dfrac{{\rm d}\omega _y }{{\rm d}t} = \dfrac{I_z-I_x }{I_y }\omega _x \omega _z + \dfrac{M_{y1}}{I_y } \\ \dfrac{{\rm d}\omega _z }{{\rm d}t} = \dfrac{I_x-I_y }{I_z }\omega _x \omega _y + \dfrac{M_{z1} }{I_z} \\ \end{array}\right\}$

(2)

对当前研究的问题：偏航方向是固定的，俯仰方向为强迫振荡，滚转方向为自由运动.于是，从动力学方程式(2)可以看到，滚转方向与俯仰方向不存在动力学方程上的耦合关系，这使我们可以专注于研究俯仰振荡过程导致的非定常流场对滚转动态特性的影响.

2 方形截面弹静态滚转气动力特性

在分析方形截面弹静态滚转气动特性之前，首先对发展的算法进行验证.这里采用航天飞机上脱落泡沫碎片的轨迹预测^[31]（图 2）验证本文发展的松耦合算法；采用方形截面弹^[32]定常计算结果验证计算程序对研究外形的模拟能力, 其中压力系数用$C_{\rm p}$表示.

图 4是不同攻角时静态滚转力矩(${ C}_1$)随滚转角变化的曲线，这里，定义"$\square$"形截面为0$^\circ$滚转角，"◇"形截面为45$^\circ$滚转角，同时定义飞行器顺时针滚转时(从底部往前看)滚转角为正.对"$\square$"形姿态，在0$^\circ\sim$12$^\circ$攻角范围内，飞行器在滚转方向是静不稳定的；攻角为10$^\circ$时，静不稳定的值最大.而在14$^\circ \sim$40$^\circ$攻角范围内，同样对"$\square$"形姿态，飞行器转变为静稳定的，且随攻角增大，静稳定性增强.静稳定和静不稳定的临界值为约13$^\circ$攻角.

对飞行器静态滚转气动特性的分析，可以通过数值模拟飞行器自由振动的过程加以验证，图 5给出了不同攻角时模拟得到的滚转角速度对滚转角的相图.飞行器从5$^\circ$滚转角位置释放后，飞行器在自身滚转力矩作用下自由振动.攻角为10$^\circ$时，飞行器的滚转运动沿正向发散，形成快速滚转运动；攻角大于14$^\circ$后，滚转运动均是振荡收敛的.

3 俯仰振荡对滚转特性的影响分析 3.1 俯仰振荡对滚转静、动态稳定性影响

首先定义强迫俯仰振荡过程为

$\alpha (t) = \alpha _0 + \alpha _{\rm m} \sin (2\pi ft)$

(3)

其中，{$\alpha $}$_{0}$为起始攻角；{$\alpha $}$_{\rm m}$为振幅，均取1$^\circ$；$f$为有量纲振荡频率，单位为Hz，当$f=0$时，表示无俯仰振荡过程.

在开始研究俯仰振荡对滚转运动特性的影响之前，首先分析其对滚转静态特性和滚转动态特性的影响.此时，保持滚转角$\gamma=5^\circ$不变，由于$\gamma=5^\circ$相当于在"$\square$"形截面位置附近偏离了一个小量，该位置的滚转力矩可以在一定程度上反映飞行器滚转方向的静态稳定性. 图 6分别给出起始攻角$\alpha $}$_{0}$为10$^\circ$，14$^\circ$，20$^\circ$和30$^\circ$时，俯仰振荡过程中滚转力矩随时间的变化曲线.

在滚转静不稳定时($\alpha $}$_{0}$=10$^\circ$)，俯仰振荡过程中的滚转力矩整体位于无俯仰振荡值的下方，降低静不稳定的量值；滚转静稳定时($\alpha $}$_{0}$=20$^\circ$和30$^\circ$)，俯仰振荡过程中的滚转力矩整体位于无俯仰振荡值的上方，降低静稳定的量值，且随着攻角增大，滚转静稳定性增强，降低静稳定的效果减弱；在临界位置($\alpha $}$_{0}$=14$^\circ$)，俯仰振荡时，滚转力矩随攻角变化围绕无俯仰振荡值上下振动，波动的幅值几乎相同.

俯仰振荡过程对滚转动态稳定性的影响可通过分析滚转力矩对攻角的相图得到. 图 7分别给出起始攻角$\alpha $}$_{0}$为10$^\circ$，14$^\circ$，20$^\circ$和30$^\circ$时，俯仰振荡过程中滚转力矩对攻角的相图.随着攻角的变化，滚转力矩对攻角的迟滞圈形状也不断改变，有正常的迟滞圈形状($\alpha $}$_{0}$=30$^\circ$)，有"8"字形($\alpha $}$_{0}$=10$^\circ$)，也有双"8"字形($\alpha $}$_{0}$=14$^\circ$).但综合来看，迟滞圈的旋转方向逆时针占优，表明俯仰振荡过程中滚转阻尼为正，飞行器对气体做功，也即会增强滚转方向的动态稳定性.

3.2 俯仰振荡对滚转运动特性的影响

在分析了俯仰振荡对滚转静、动态稳定性的影响之后，再放开滚转方向的自由度，研究俯仰振荡对滚转运动的影响.俯仰振荡的形式仍按照式(3)给定，图 8给出了不同起始攻角时，俯仰振荡对滚转运动影响的模拟结果.由于飞行器发生连续快速滚转时，滚转角的变化范围较大，为便于分析，在绘图时通常将滚转角进行转换.滚转角变化范围定为$\gamma $=0$^\circ$至$\gamma $=-360$^\circ$，超出范围后将滚转角转换到该范围中.因此图 8中出现一条虚线段连接两条实线段的端点，即表明飞行器绕体轴滚转了一周，即360$^\circ$.

对图 8(a)，起始攻角$\alpha _{0}$=10$^\circ$.在10$^\circ$攻角时，飞行器的滚转方向是静不稳定的，释放后飞行器将发展为连续滚转运动模式（参考图 5）.但在俯仰振荡的影响下，飞行器的滚转运动规律有很大不同.当俯仰振动的频率较低时（$f=5$ Hz），飞行器的滚转运动仍会发展为快速滚转运动；但由于俯仰振荡过程会降低静不稳定度，增强动稳定性，振荡越快，这种效果越明显，因此随着俯仰振荡的频率升高，滚转运动逐步转化为震荡收敛形式，最终会稳定在"$\square$"形截面位置.

对图 8(c)和图 8(d)，起始攻角$\alpha_0$分别为20$^\circ$和30$^\circ$.此时飞行器滚转方向是静稳定的，强迫俯仰振荡过程会降低静稳定的值，但会增强动稳定性.从模拟结果来看，$\alpha _{0}$=20$^\circ$时，滚转的静稳定性相对较低，强迫俯仰振荡过程对滚转角收敛的位置有较大影响；随着$\alpha _{0}$增大到30$^\circ$，滚转的静稳定性增强，这种影响已经很小.

而对图 8(b)，起始攻角$\alpha _{0}$=14$^\circ$，情况则复杂很多.飞行器本身处于静稳定和静不稳定的临界点，不同俯仰振荡频率的模拟结果几乎没有规律可循：滚转角既有可能收敛，不同俯仰振荡频率时收敛的位置也不同；也有可能沿正向或负方向发散；甚至也有滚转两周后又转变为收敛的运动形态.

这些结果也表明，俯仰振荡过程对滚转运动特性（收敛状态，收敛位置等）的影响较为复杂，具体影响程度视滚转静稳定情况而定.

3.3 网格无关性研究

在以上分析中，虽然通过定常计算结果对程序进行了验证，但由于计算网格较稀，仍可能对非定常的计算结果产生较大影响.为此，需开展网格无关性研究.

这里采用三套计算网格G1，G2和G3，计算网格的拓扑结构相同，均为多块对接网格，网格点数量分别约为450万、880万和1 400万. 图 9给出了G1和G3网格在对称面上的网格分布比较，相对于G1网格，G3网格在流向、法向和周向的网格点数均有增加.

俯仰振荡对滚转运动影响的网格无关性研究选取初始攻角$\alpha _{0}$=10$^\circ$，振荡频率为20 Hz的计算状态，计算条件同上. 图 10给出了3套网格的计算结果比较，可以看到：首先，随着网格数量的增加，计算结果的差异性不断降低；其次，网格确实对计算结果有较大影响，但不同网格计算得到的俯仰振荡对滚转运动的影响规律基本一致.

4 结论

基于建立的六自由度刚体动力学方程和N-S方程的耦合模拟技术，用数值方法研究了方形截面弹俯仰振荡过程对滚转特性的影响，主要结论有：

(1) 研究的方形截面弹外形，其外形特点决定了滚转方向与俯仰方向不存在动力学方程上的耦合关系，这使得可以专注于研究俯仰振荡过程导致的气动力非定常效应对滚转特性的影响，这与传统上研究的三角翼等面对称外形的俯仰/滚转耦合特性有很大不同.

(2) 定性分析了俯仰振荡过程对滚转静态稳定性和动态稳定性的影响.俯仰振荡过程一般会增强滚转动态稳定性；滚转方向为静稳定时，俯仰振荡过程会降低静稳定度；滚转方向为静不稳定时，俯仰振荡过程则会降低静不稳定度.

(3) 分析了俯仰振荡过程对滚转运动特性的影响.对滚转静不稳定情况，快速的俯仰振荡也可能使飞行器滚转运动呈收敛状态；对滚转静稳定情况，静稳定性较弱时，俯仰振荡对滚转运动收敛的位置有较大影响，随着静稳定性增强，这种影响也随之减弱.