窄条翼导弹模型摇滚运动动力学特性研究

达兴亚，周为群，赵忠良，陶 洋

（中国空气动力研究与发展中心 高速空气动力研究所，绵阳 622661）

0 引 言

高机动性导弹突出强调大迎角快速转弯能力，极小展弦比窄条翼在大迎角不存在涡破裂现象，纵向气动特性极为优异，是高机动性导弹经常采用的布局形式［1］（图1）。但是，窄条翼导弹强烈的非线性会带来严重的翼／舵气动干扰，导致全弹气动特性随滚转角和迎角急剧变化，在机动时极易诱发非指令的自激摇滚振荡。

图1 窄条翼示意图（NASA和德国IRIS－T）Fig.1 Strake missiles（NASA and German IRIS－T）

摇滚是飞行器滚转方向的自由振荡，由于以滚转方向为主，一般简化为单自由度运动来开展研究［2］。通常有三种途径：一是通过风洞实验获取摇滚特性；二是通过数值模拟摇滚运动，分析流场变化，从微观上搞清摇滚的流动机制；三是通过非线性动力学理论，分析摇滚运动的动力学特性，搞清摇滚的产生机理，给出失稳判据和预测方法。张涵信等人指出［3］，搞清摇滚的产生机理是十分必要的：即什么条件下摇滚是稳定的，什么条件下演化为周期摇滚。由于国内研究起步晚，对这类导弹的摇滚机理还了解不多。

在过去的工作中［4］，对窄条翼布局导弹模型进行了风洞摇滚实验，发现窄条翼导弹模型具有两个摇滚区间，第一摇滚区间在20°迎角附近，迎角范围极小；第二摇滚区间在35°迎角，迎角范围较宽。实验同时明确了摇滚运动与窄条翼／尾舵干扰密切相关，但尚不清楚对应的动力学特性。本文针对模型第二摇滚区开展数值计算，重点分析摇滚的动力学特性，搞清摇滚的产生机理。

1 数值计算方法

1.1 控制方程

摇滚数值模拟常使用的控制方程有速度势方程［5］、欧拉方程［6］，以及 NS方程［7］。对三角翼研究表明，速度势方程和欧拉方程能成功模拟出极限环振荡，同时也表明摇滚可能主要是非定常无粘现象主导［8］。三维NS方程对分离预测更准确，能得到与实验更接近的结果。本文流场控制方程采用三维非定常NS方程：

式中Q为守恒变量，F、G、H为无粘矢通量，Fv、Gv、Hv为粘性矢通量。转动方程为定义在弹体坐标系的刚体运动方程：

其中φ为滚转角，Ix为转动惯量，ω为角速度，M为力矩，上标b表示弹体坐标系下的分量。

1.2 计算方法［9］

流场求解器是基于结构网格的有限体积求解器，粘性项采用James中心差分，无黏项采用Roe格式，以F为例：

各参数具体意义参见文献［10］。时间推进采用双时间步法［11］，湍流模型选取S－A 模型［12］。

转动方程使用Adams预估校正法求解，气动／运动耦合计算采用双时间步Adams预估校正策略，具体思路是：在CFD完成n时刻计算并得到气动力后，采用Adams－Bashforth预测格式预测n＋1时刻的滚转角和滚转角速度，更新姿态和CFD边界条件；之后，CFD采用双时间步法开始计算n＋1时刻流场，在双时间步的每一内迭代步结束后，CFD得到新的气动力，此时采用Adams－Moulton格式修正滚转角和角速度，在内迭代结束或者状态变量收敛后，推进到n＋2时刻。

在双时间步内迭代过程实现隐式修正，三阶Adams－Moulton格式为：这里n代表时间步，N代表内迭代步，如果计算收敛，随着内迭代步的推进，xn＋1，N＋1将趋于常值。

研究表明，以上三阶耦合计算方法在保证流场收敛和一定计算精度条件下，可以显著增大时间步长，缩短计算时间，详细分析参见文献［9］。

2 摇滚运动数值模拟结果及讨论

2.1 计算网格及状态

计算模型为钝头／窄条翼／尾舵组合导弹模型（类似IRIS－T，如图1），网格为标准多块结构对接形式，在模型的周向分布了223个网格点，尾舵网格进行了加密处理，网格总量约580万，其上游距头部10L，下游距后缘10L，远场边界距中心线7L，L为全弹长度。为提高并行计算效率，将网格分为111块，最大块网格量约15万，保证每个计算核心分配到大致相等的计算量。

计算针对窄条翼导弹模型的第二摇滚区，Ma＝0.6，α＝35°，“×”字型状态，Re和Ix取风洞试验对应值。计算采取预偏滚转角φ0、初始角速度为0的起动策略。

图2 计算网格空间拓扑Fig.2 Space topo of computational grid

2.2 计算结果及讨论

由图3可见，模型从φ0＝5°自由释放，在经历短暂的几个振荡周期后，模型很快进入极限环振荡（Limit Cycle Oscillation，LCO），振幅约10.14°，频率20Hz。图4力矩迟滞曲线形状上比较“肥胖”，说明迟滞现象明显，模型吸收和耗散能量的速度很快，导致模型很快进入极限环状态。图5相位图可以看出，进入极限环振荡后，新周期的相位与前一周期完全重合。

图3 滚转角时间历程Fig.3 Time history of roll angle

图4 迟滞曲线Fig.4 Hysteresis loop

图5 相平面Fig.5 Phase plane

此状态下对应的风洞试验结果为：振动均方差σγ＝8°，频率14.3Hz，计算结果与试验结果振幅吻合较好，而频率比试验值大5.7Hz。引起频率差异主要有两方面原因：一是在气动载荷的作用下，模型转动时可能存在较大的机械阻尼，目前在计算中还没有一种普适的方法引入机械阻尼；二是模型转动时连轴承一起转动，计算时没有计入这部分转动惯量。

极限环显著的特征是一个周期内吸收与耗散的能量相等，总能量变化为0［13］。当模型处于低能状态时，会不断吸收能量而趋向极限环，当处于高能状态时，会不断耗散能量而趋向极限环。图5中由内向外发展的相位图即是由低能态向极限环过渡的过程，图4中的内环为顺时针方向，系统吸收能量，两边环是逆时针方向，能量被耗散，其面积和等于顺时针环面积，将吸收的能量耗散掉，一个周期内系统能量变化为0，振荡既不发散，也不收敛，维持极限环状态。

3 摇滚动力学特性分析

3.1 摇滚产生的动力学机理

图6给出了迎风面和背风面窄条翼和尾舵的迟滞曲线。从斜率上看，背风舵和背风翼的斜率为负，且斜率较大，迎风舵迟滞环虽大，但斜率小，迎风翼的斜率为正，因此，导弹的横向静稳定性主要由背风舵和背风翼提供。

四条曲线中仅迎风舵力矩为顺时针环，即有能量注入，运动发散；背风翼和背风舵为逆时针环，系统对外做功，能量被耗散。迎风翼也为逆时针单环，虽然面积较小，但对系统的动稳定性贡献为正。也可以看出，迎风舵在平衡滚转角处的迟滞环面积很大，而背风舵的迟滞环面积很小，导致背风舵不能抵消迎风舵的动不稳定性，使模型丧失滚转阻尼，吸收能量，运动逐渐发散；随着振幅的增大，背风舵的动稳定性增强，而迎风舵的动不稳定性降低，整个模型从动不稳定变为动稳定，耗散能量，运动受到抑制，最终进行等幅等周期极限环振荡。

图6 窄条翼和尾舵迟滞曲线Fig.6 Hysteresis loops for strakes and fins

3.2 初始滚转角对摇滚的影响

大多数情况下，极限环是稳定的，模型从任意非平衡状态都将稳定到该极限环［14］；特别情况下极限环是不稳定的，系统还存在其它极限环，模型在扰动作用下会在极限环之间来回跳动，进入混沌运动［15］。

取初始滚转角γ0＝15°，模型同样进入极限环振荡，如图7所示，相图由外向内发展，即由高能态过渡到极限环。该极限环与γ0＝5°极限环基本重合，说明确实存在一个稳定的极限环，模型在高／低能状态或初始微扰动作用下，都将进入该极限环振荡。在模拟时未考虑外部扰动，如果外部扰动过大，极限环的稳定性不再起作用，也可能进入混沌运动，这种情况不在本文讨论范围。

图7 初始滚转角的影响Fig.7 Effect of initial roll angle

3.3 转动惯量对摇滚的影响

摇滚实验采用了铝制模型，未作转动惯量影响研究，为了明确转动惯量的影响，假设采用钢制模型（密度约为铝的3倍），将转动惯量增大至3倍后进行模拟。模拟的振幅9.3°，频率9.4Hz，振幅变化不大，频率变化明显。

由图8中力矩迟滞曲线可见，力矩迟滞曲线都呈双“8”环，环的大小基本一致，只是两边的小环稍有差异。比较图5和图9，增大转动惯量导致最大滚转角速度从1200°／s降低到600°／s，说明在这一较宽广的转速范围，气动和运动的耦合特性基本相似，角速度对力迟滞特性和流场跟随性的影响较小。

图8 转动惯量的影响Fig.8 Effect of moment of inertia

图9 相平面（Ix＝3I0）Fig.9 Phase plane（Ix＝3I0）

4 结 论

通过对窄条翼导弹模型的摇滚运动模拟，得到了摇滚的动力学特性，结论如下：

1）多个模拟状态表明，窄条翼导弹模型力矩迟滞曲线呈双“8”环，中间环为顺时针方向，两边环为逆时针方向；

2）窄条翼导弹模型摇滚的动力学原因是：由于背风舵不能提供足够的动稳定性，模型在迎风舵的动不稳定作用下运动发散，随着振幅的增大，背风舵动稳定性增强，迎风舵动不稳定性减弱，从而抑制了振荡发散，最终进入极限环振荡；

3）第二摇滚区的极限环是稳定的，模型在任意初始状态或微扰动作用下都将进入该极限环振荡；

4）在非定常特性较强时，存在一个转动惯量范围，摇滚振幅变化不大，振荡周期变化明显。

［1］ALLEN J M.Aerodynamics of an axisymmetric missile concept having cruciform strakes and in－line tail fins from mach 0.60to 4.63［R］.NASA TM 2005－213541，2005.

［2］KATZ J.Wing／vortex interactions and wing rock［J］.Progress inAerospaceSciences，1999，35：727－750.

［3］ZHANG H X，LIU W，XIE Y F，et al.On the rocking motion and its dynamic evolution of a swept delta wing［J］.ACTA AerodynamicaSinica，2006，24（1）：5－9.（in Chinese）张涵信，刘伟，谢昱飞，等.后掠三角翼的摇滚及其动态演化问题［J］.空气动力学学报，2006，24（1）：5－9.

［4］DA X Y，ZHAO Z L，TAO Y，et al.Experimental investigation on free－to－roll motion of strake missile experimental［J］.JournalofExperimentsinFluidMechanics，2012，26（6）：40－44.（in Chinese）达兴亚，赵忠良，陶洋，等.窄条翼导弹模型摇滚特性试验研究［J］.实验流体力学，2012，26（6）：40－44.

［5］KONSTADINOPOULOS P，MOOK D T，NAYFEH A H.Subsonic wing rock of slender delta wings［J］.JournalofAircraft，1985，22（3）：223－228.

［6］KANDIL O A，SALMAN A A.Three－dimensionla simulation of slender delta wing rock and divergence［R］.AIAA－92－0280，1992.

［7］LIU W，ZHANG H，ZHAO H.Numerical simulation and physical characteristics analysis for slender wing rock［J］.JournalofAircraft，2006，43（3）：858－861.

［8］ARENA A，NELSON R.A discrete vortex model for prediction wing rock of slender wings［R］.AIAA－92－4497，1992

［9］DA X Y，TAO Y，ZHAO Z L.Research on virtual flight simulation based on prediction－correction coupling method and chimera grid［J］.ActaAeronauticaetAstronauticaSinica，2012，33（6）：977－983 （in Chinese）达兴亚，陶洋，赵忠良.基于预估校正和嵌套网格的虚拟飞行数值模拟［J］.航空学报，2012，33（6）：977－983.

［10］ROE P.Approximate Riemann solvers，parameter vectors，and difference schemes［J］.JournalofComputationalPhysics，1981，43（2）：357－372.

［11］DUBUC L，et al.Solution of the unsteady Euler equations using an implicit dual－time method［J］.AIAAJournal，1998，36（8）：1417－1424.

［12］SPALART P R，ALLMARAS S.A one－equation turbulence model for aerodynamic flows［R］.AIAA－92－0439，1992.

［13］NELSON R C，PELLETIER A.The unsteady aerodynamics of slender wings and aircraft undergoing large amplitude maneuvers［J］.ProgressinAerospaceSciences，2003，39：185－248.

［14］LIU Z L.Control and simulation investigation of wing rock phenomenon［D］.［PhD thesis］.Concordia University，Canada，2004.

［15］RONG Z，DENG X Y，WANG B，et al.Experimental investigation on flow characteristics during double limit cycle oscillation induced by forebody vortex［J］.JournalofExperimentsinFluid Mechanics，2010，24（3）：19－24.（in Chinese）荣臻，邓学蓥，王兵，等.前体涡诱导双极限环摇滚流动特性的实验研究［J］.实验流体力学，2010，24（3）：19－24.