正规形法在弹箭非线性运动分析中的应用*

李东阳，常思江，王中原，魏 伟

(1. 南京理工大学 能源与动力工程学院， 江苏 南京 210094； 2. 瞬态冲击技术重点实验室， 北京 102202)

弹箭非线性运动理论始于20世纪50年代，由于其对工程应用的巨大指导作用，一直以来都是外弹道理论的研究重点[1-3]。已有研究表明，在某些特定条件下，非线性效应对弹箭运动特性的影响不容忽视，甚至起决定性作用[4-6]。某些基于线性理论设计的弹箭在飞行试验中，会发生预料之外的质心运动或姿态变化，无法完成预定的飞行任务。例如：无控十字翼导弹的转速在一定条件下被锁定在某一转速附近而无法达到设计值，同时，攻角也被锁定而产生圆锥运动，此种现象被称为转速-攻角闭锁[4, 7]。此时，如果被锁定的转速恰好位于该弹箭的共振转速区域，则可能飞行失稳甚至掉弹。对于有控弹箭，攻角和转速被锁定将大大增加飞行控制的难度，使其往往无法达到良好的控制效果[6，8-10]。

研究表明，引起弹箭非线性运动的因素有很多，如气动力非线性、结构非线性、几何非线性等，其中气动力非线性占主导作用[11-12]。有学者发现[5]，一些不常考虑的气动力和力矩往往可用于解释某些非线性现象。也有不少学者致力于提高气动力和力矩模型的精确性(描述成攻角的高阶多项式)，以期更接近实际地描述弹箭非线性运动。Liano等[13]研究发现，为了准确预测转速闭锁现象，在攻角为12°和20°时，与滚转角有关的非线性力矩关于攻角的阶数应分别不得低于五次和七次；Morote[9]研究发现，对于某种弹箭结构，当静力矩系数关于攻角的阶数达到七次时才能准确地拟合试验数据。显然，为了更好地研究弹箭非线性运动，引入高阶气动力系数是十分必要的。

当前，制约弹箭非线性运动理论进一步发展的主要因素是缺乏有效的方法，当攻角方程中存在高阶气动力系数及其耦合项时，获得高精度的解析解相当困难。传统的弹箭非线性研究方法[1]都存在一定的局限性，如：奇异摄动法对弱非线性系统原则上能给出满足任意精度要求的解析解，但高阶精度计算量大，解的形式过于繁杂，难以用于运动特性的进一步分析；而简化计算的周期平均法作为一次近似理论，只限于定性研究，不能用于对精度要求高的定量计算，因此适用于弱非线性系统；渐近法虽然在一定程度上融合了上述两种方法的优点，但随着精度的提高，复杂度迅速增加。在弹箭非线性运动学研究中，角运动通常运用拟线性法进行分析，即假设在非线性情况下，弹箭的角运动仍能写成线性运动解的形式，再利用周期平均法求出可变的频率和阻尼指数，并利用振幅平面方程和奇点理论求得稳定性条件。缺点是静力矩非线性较强(如考虑三次非线性静力矩)时误差较大。虽然在考虑三次静力矩时能借用椭圆函数得到精确解，但椭圆函数无法处理三次以上的静力矩项。外弹道学领域中常用的摄动法，是在以三次静力矩作用下椭圆函数表示的精确解的基础上，将其他非线性力矩的影响作为椭圆函数基础解的小扰动，但也仅限于弱非线性的处理[1]。以上方法在应用上的局限性已为人所共知，因此不断有研究探讨新的攻角方程解法和稳定性分析方法。文献[14]尝试用同伦分析法求解高次非线性静力矩作用下的角运动方程，得到了有效的解析解和相对保守的稳定攻角范围，但为了确保解的收敛需引入特定算法。霍普夫分岔理论也被用于分析弹箭非线性运动的稳定性[15]。文献[16]利用推广的打靶法计算了角运动周期解的幅值和周期,结合Floquet理论分析了周期解的稳定性。

正规形法的基本思想是通过引入恒等变换，将原系统微分方程化为尽可能简单的形式，它不仅能够在更大区间内给出非线性方程的有效解析解，同时能较容易地获得非线性系统的频率随幅值以及振动随自变量的变化曲线[17-19]。频率随幅值的变化曲线在振动模态分析中发挥了重要作用。而幅值变化曲线，则有助于分析运动稳定性，这是前述奇异摄动法、平均法以及多尺度法难以给出的结果。利用正规形法得到系统的正规形，可用于进一步分析系统在参数摄动情况下的分叉现象，但暂不在本文讨论范围内。近来，有许多对正规形法本身的研究和拓展[18, 20-23]，同时其在电力系统分析[20, 24-25]、振动模态分析[18, 23, 26]、飞行器稳定性分析[27]、生态研究[28]、转子动力学[29-31]等领域均得到了广泛应用，成为研究非线性微分方程的有力方法。而从现有文献来看，正规形法在外弹道学领域中鲜有应用。

因此，本文的主要目的和出发点就是对弹箭非线性运动分析问题引入新的解析方法，探讨正规形法对弹箭非线性运动研究的优势。本文首先简要介绍正规形法的构造原理，在此基础上，将针对具有高阶静力矩系数(五次和七次)的弹箭攻角运动方程，在有阻尼(二次)和无阻尼条件下，利用正规形法解析求解并进行稳定性分析，用数值积分验证结果的可行性，以期为弹箭非线性运动理论研究提供新的有效方法。

1 正规形构造方法

考虑非线性系统

(1)

式中：自变量x=[x1,x2,…,xn]T；A为n×n常数矩阵；n×1向量F(x)为非线性项，在原点处至少二阶连续可微。为不失一般性，假设其平衡点为原点或已通过线性变换将感兴趣的平衡点移至原点。

正规形法的基本思想是采用一定的近似恒等坐标变换，将复杂的动力学系统转化为其等价类中形式最简单的那一个，以便进一步分析和求解。简化通常在系统动态平衡点或周期轨道(如极限圆)的邻域内进行[17]。由于所考虑的弹箭角运动方程的形式特殊，A具有一对纯虚特征根，这里仅简单介绍所用的特殊构造方法，即复数表示法[17]。

引入线性坐标变换x=Pz得到式(1)的Jordan形式为

(2)

(3)

引入复变量

(4)

(5)

其中，i为虚数单位。这样，与正规形基本构造方法相比，在得到相同正规形形式的同时省去了许多复杂的代数运算。此时，近似恒等变换可选为

(6)

共振条件为

m1-m2=1

(7)

2 攻角平面内弹箭非线性运动的解析解

求出弹箭角运动方程的解析解，对弹箭运动和稳定性分析都十分方便，尤其是当前各种新型弹箭层出不穷，气动特性和运动特性更为复杂，高阶气动力和力矩的研究变得更加重要。已有研究表明[13]，传统的三次形式已不足以准确刻画静力矩，试验数据表明，即使在15°的中等攻角，取至攻角的七次也不足以准确刻画静力矩中的高度非线性。高次非线性项对弹箭所受静力矩的准确刻画非常重要，而静力矩作为重要的气动力矩，研究其对弹箭角运动的影响就变得十分重要。此外，当二次非线性阻尼力矩同时存在时，会诱发圆锥运动，在弹箭平面运动中显示为平面等幅振荡[1-3, 32]。而当非线性力矩较强时，传统拟线性法与平均法的结合将无法得到满足精度的结果。因此，本文将在传统仅考虑三次非线性静力矩的基础上，采用正规形方法研究五次和七次非线性静力矩以及二次非线性阻尼对角运动的影响。

弹箭的平面角运动方程可描述为

δ″+(H0+H2·δ2)δ′-(M0+M2·δ2+
M4·δ4+M6·δ6)δ=0

(8)

不难看出，δ=0、δ′=0为该平面角运动系统的一个平衡点。注意，对于不同的气动力系数组合，也可能存在其他有意义的实平衡点。在实际应用中，研究人员往往对弹箭在原点周围的角运动更感兴趣。因此，利用正规形理论寻求七次非线性角运动方程在原点附近的正规形式，并研究不同气动系数组合对攻角幅值、振动周期以及稳定性的影响。

2.1 考虑高阶静力矩的方程解析解

(9)

式中：f(δ,δ′)=-(H0+H2·δ2)δ′+(M0+M2·δ2+M4·δ4+M6·δ6)δ

记攻角为

利用式(5)～(6)，为求简洁，详细推导不再赘述，可得当角运动方程取得最简形式时，近似恒等变换式(6)的系数应选择如下

(10)

由此，可保证消去所有非共振项后得到尽可能简单的正规形式。式中的“？”项表示未在方程展开式中出现的待定系数Γli，故无法确定此系数的具体值，即可为任意值。为得到最简单的正规形，将其取为0。此时，可得角运动方程的正规形为

(11)

通过逆变换可得出攻角的具体表达式。首先，给出正规形式(11)的平凡解，即

(12)

δ=acos(Ωs+β0)+

(13)

其中

(14)

Ω=ω0+β′

(15)

其中：k=1,3,5,7,…；n=3,5,7,…且n≥k；k1=1,3,…；n1=1,3,…。

分析该表达式，可知：①攻角中的余弦部分系数只和静力矩系数有关，而和阻尼系数无关；而正弦部分系数则只和阻尼系数、线性静力矩系数M0有关，而和高阶非线性静力矩无关。②正弦部分的幅值远小于相应余弦部分，且包含阶数较低的谐波。

为全面确定攻角表达式，将平凡解式(12)代入正规形式(11)，展开后将实部和虚部分开，得

(16)

(17)

也适用于七次以上的非线性静力矩。

由式(16)可知，攻角幅值线性部分的变化特性由阻尼系数H0、H2决定。当H0=0、H2=0时，a′=0，即无阻尼时，弹轴若稳定摆动则幅值恒定为初始角运动对应的幅值。但要注意，该正规形是在原点附近推导的，当角运动方程式(8)存在其他实平衡点时，这些平衡点的稳定性将影响弹轴维持稳定等幅摆动的初始攻角范围。

对于弹道弧长s积分，可得线性幅值a(s)随s的变化规律为

(18)

式中：a0由初始条件δ(0)=δ0,δ′(0)=δ′0确定

(19)

值得注意的是，该式也可用于确定稳定初始条件范围，但由于该式为线性近似关系，实际使用时为保证正确性可做一定的缩放。此外，若将式(18)代入式(17)，可得弹轴摆动频率随弹道弧长s的变化规律Ω(s)。因此，攻角δ可完全表示为弹道弧长s和初始条件δ0、δ′0的函数。

式(17)体现的是攻角频率受幅值影响的非线性特性。频率主要受静力矩系数影响而和阻尼项无关，结合式(17)可分析不同静力矩系数对弹轴摆动频率的影响，如Mi<0，弹轴摆动频率必然增大；而Mi>0，弹轴摆动频率必然减小。

2.2 稳定性分析

2.2.1 无阻尼时角运动稳定性分析

由于系统的平衡点δ*满足

(20)

因此，在高阶非线性静力矩系数影响下，角运动可能存在除原点以外的其他0～3对实平衡点。平衡点对应的特征根λ满足

λ2=M0+3M2δ*2+5M4δ*4+7M6δ*6

(21)

故原点对应的特征根满足λ2=M0<0，据线性化原理，其为中心。

2.2.2 阻尼对角运动的影响

阻尼项的存在虽不改变角运动方程(8)的平衡点位置，但影响其稳定性，此时特征根满足

λ2+(H0+H2δ*2)λ-(M0+3M2δ*2+

5M4δ*4+7M6δ*6)=0

(22)

为了确定攻角运动的稳定性，需综合考虑系统平衡点式(20)和幅值方程(16)。原点的稳定性将由线性静力矩系数M0和线性阻尼系数H0共同确定，当M0<0时，若H0<0则原点为不稳定平衡点；若H0>0则原点为稳定平衡点。

考虑幅值方程(16)，H0、H2不同符号组合决定了不同的幅值变化规律，如表1所示，可能存在稳定或不稳定极限环。

表1 幅值方程(16)平衡点分析

表1所示情形1条件下，原点变为不稳定平衡点，且可能存在极限环δLC=-4H0/H2，但并不能保证弹轴最终以幅值δLC进行等幅摆动。实际上，在某些气动系数组合下弹轴也可能完全失稳。以三次非线性静力矩为例进行分析。

对于表1所示情形3，原点为不稳定平衡点，若M2<0，则不存在其他平衡点，任意初始条件下，攻角都将失稳；若M2>0，由特征根方程(22)知，平衡点均为不稳定平衡点，因此攻角也都将发散。

综上可见，高阶静力矩系数导致角运动方程的平衡点个数和稳定性改变，考虑二次阻尼项时如果极限环存在，则角运动的稳定性将由两者共同确定。因此，仅用传统的平衡点分析法无法确定非线性角运动的具体运动规律和稳定初始条件，而结合正规形分析给出的幅值方程(16)则给出了较为准确的分析。

3 解析解的精度验证及相关分析

本节将通过数值积分验证正规形给出的角运动解析解及稳定性分析的正确性，为便于对比，这里采用文献[14]中给出的气动参数组合，如表2所示。

表2 数值积分验证用气动参数[14]

3.1 考虑无阻尼五次静力矩的解析解

3.1.1 稳定初始条件范围的确定

无阻尼情况下，首先计算四种气动力条件下的攻角平衡点δ*及其对应特征根，结果如表3所示。

表3 无阻尼五次静力矩角运动方程稳定性分析

由表3结果可知，情形Ⅰ和情形Ⅱ均为全局稳定；情形Ⅲ和情形Ⅳ均存在原点之外的平衡点且为鞍点，由于原点为中心，因此相平面轨迹分别如图1(a)、图1(b)所示，取k=0.90、C(0.90δ*)与数值积分确定的最大初始条件范围CM相比，正规形方法给出的初始条件关系式(19)，虽相对保守，但可较快捷地确定稳定初始条件范围。

(a) 情形Ⅲ(a) Case Ⅲ

3.1.2 正规形攻角解有效性验证

下面选取不同的初始攻角δ0(表3中的情形Ⅰ～Ⅳ分别对应8°、 10°、 19°、 6°)，在初始攻角速度δ′0=0时，分别用正规形法给出的攻角解析解和数值积分方法计算攻角随弹道弧长的变化，如图2所示。在一定的初始攻角范围内，正规形解析解和数值积分结果(下同)吻合得很好，这表明采用正规形法获取五次静力矩作用下攻角运动方程解析解的有效性。

(a) 情形Ⅰ(a) Case Ⅰ

图3给出了四种气动组合下，弹箭稳定摆动周期T=2π/Ω与攻角线性幅值a之间的关系。其中，线型对应正规形解析解，符号标志对应数值解。当攻角不大时(如小于15°)，正规形解析解较好地预测了频率与幅值之间的关系；当攻角较大时，情形Ⅲ的结果依然吻合得较好，其余三种情况存在一定误差，但总体上仍在可接受的范围内。这主要是由于a仅为攻角幅值的线性部分，当攻角较大时，非线性部分的影响变强，故用a预测出的周期会与实际值产生一定的偏差。如图4所示，由实际攻角初值δ0与a0之间的关系可见，对于情形Ⅲ，攻角大于20°时，δ0与a0相差较大，因此攻角小于20°时，直接用δ0代替a0得到的解析解与真实解吻合得较好。

图3 周期和线性幅值关系Fig.3 Period variation with amplitude

图4 正规形法得出初始攻角δ0与初始线性幅值a0之间的关系(情形Ⅲ)Fig.4 Relationship between δ0 and a0 by method of normal forms (CaseⅢ)

3.2 考虑有阻尼五次方静力矩的解析解

3.2.1 极限环的形成及对稳定性的影响

仍采用文献[14]算例中的气动力系数：M0=-5×10-5,M2=-4.5×10-5,M4=-8×10-3，组合不同的阻尼系数得到表4所示的气动系数组合。其中，情形Ⅴ为无阻尼，情形Ⅴ.1为稳定极限环，情形Ⅴ.2为不稳定极限环。这些组合情形将用于验证正规形解析解在有阻尼条件下的有效性。

表4 有阻尼五次方静力矩角运动方程平衡点稳定性分析

从表4和表3可以看出，情形Ⅴ和情形Ⅰ相似，仅有原点一个平衡点，因此弹轴均以初始条件确定的幅值做稳定周期摆动，如图5和图6所示。而在有阻尼的情况下，如情形Ⅴ.1，负阻尼使得原点变得不稳定，于是当初始攻角条件在原点附近时，攻角将会发散，但正阻尼H2的存在避免了攻角无限增大，最终稳定在δLC=8.10°，该稳定极限环也吸引了环外的轨道，使得初始条件在环外的轨道脱离自由振荡而最终稳定于极限环，如图5(a)所示；图5(b)～(c)给出了情形Ⅴ.1在不同初始条件下正规形解和数值积分结果，两者吻合得很好。对应不稳定极限环的情形Ⅴ.2则恰好相反，其相平面图如图6(a)所示，由于极限环排斥环内、环外轨道，当初始条件处于极限环内时，攻角收敛到0，处于极限环外时则发散。

(a) 相平面轨迹(a) Phase trajectory

(a) 相平面轨迹(a) Phase trajectory

图6(b)～(c)给出了情形Ⅴ.2在不同初始条件下正规形解和数值积分的结果，结果表明，两者吻合得很好。综上可见，正规形解很好地预测了攻角收敛和发散的趋势，在图示攻角范围内，线性幅值仍占主导作用。

3.2.2 极限环与平衡点的共同作用

图7 情形Ⅲ.1不稳定平衡点对极限环吸引域的影响Fig.7 Influence of unstable equilibrium point on the attraction region of limit cycle in case Ⅲ.1

对比无阻尼情形Ⅲ所对应的图1(a)与此处有阻尼情形的图7可知，阻尼的存在使原点由稳定中心变为不稳定焦点，从而造成了相轨迹的变化。

对比情形Ⅴ.1全局收敛于极限环，此处情形Ⅲ.1中其他不稳定平衡点的存在限制了周期轨道的吸引域。实际上，在弹箭角运动过程中，全局稳定的可能性不大，这是因为当攻角较大时，将会产生强烈的非线性从而使气动系数产生较大变化，极限圆的幅值和平衡点的位置、稳定性都将发生改变而形成新的攻角运动形态，原周期轨道有可能不复存在。

3.3 考虑无阻尼七次静力矩的解析解

在情形Ⅰ的基础上，本文合理假设七次非线性静力矩系数取值M6=-0.01，形成如表5所示的气动参数组合，用于验证正规形法求取七次静力矩角运动方程解析解的有效性。

表5 有阻尼七次静力矩角运动方程稳定性分析

如图8和图9所示，对于无阻尼情形Ⅵ和有阻尼情形Ⅵ.1，正规形解析解和数值积分结果吻合得很好，说明正规形解对七次静力矩角运动方程依然是有效的，并准确地预测出表4情形Ⅴ.1作用下稳定极限环的存在。由于攻角稳定在8.10°，因此解析解可以在更大的攻角范围内给出较好的结果。

(a) δ0=5°

(a) δ0=5°

七次非线性静力矩下的稳定性分析和五次类似，只是其他可能存在的实平衡点从最多2对增加到3对。若这些平衡点与可能存在的极限环幅值接近，则在分析角运动稳定性时应注意考虑系统平衡点和极限环的相互作用。若形成的极限环或平衡点都较大，超出了现实中可能存在的攻角范围，则只需分析原点的稳定性就足以确定攻角运动的稳定性，此时，正规形方法给出的解析解将完全能够满足计算或分析的精度要求。

4 结论

本文将正规形法引入弹箭非线性运动理论研究领域，通过求解一个考虑阻尼项和七次静力矩项的弹箭攻角运动方程，展示了正规形的构造过程，得到了形式简洁且通用的攻角解析表达式，也适用于更高次的非线性静力矩。与数值积分对比，验证了常用攻角范围(如小于20°)内，该解析解和数值仿真吻合得很好：无阻尼情况下，攻角较大(如大于20°)时相位误差较大；而有阻尼情况下，由于收敛到零或极限环的存在，能在更大的攻角范围内吻合得很好。同时，本文利用正规形解析解，给出了弹轴幅值和摆动频率与弹道弧长之间的关系，综合考虑幅值方程与攻角解析解的初始条件，给出了较简洁的稳定初始条件范围计算方法；通过将幅值方程和角运动方程相结合，分析了阻尼项对角运动稳定性的影响，评估了极限环的产生条件和稳定范围。上述研究结果表明，正规形方法为弹箭非线性角运动特性分析提供了一个有力的方法，后续将开展更为深入的应用研究。