正規形法在彈箭非線性運動分析中的應用*

李東陽，常思江，王中原，魏 偉

(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°)時相位誤差較大；而有阻尼情況下，由于收斂到零或極限環的存在，能在更大的攻角范圍內吻合得很好。同時，本文利用正規形解析解，給出了彈軸幅值和擺動頻率與彈道弧長之間的關系，綜合考慮幅值方程與攻角解析解的初始條件，給出了較簡潔的穩定初始條件范圍計算方法；通過將幅值方程和角運動方程相結合，分析了阻尼項對角運動穩定性的影響，評估了極限環的產生條件和穩定范圍。上述研究結果表明，正規形方法為彈箭非線性角運動特性分析提供了一個有力的方法，后續將開展更為深入的應用研究。