鐘揚威，王良明，傅健，常思江

（南京理工大學能源與動力工程學院，江蘇南京210094）

彈箭非線性角運動穩(wěn)定性Hopf分岔分析

鐘揚威，王良明，傅健，常思江

（南京理工大學能源與動力工程學院，江蘇南京210094）

為了分析彈箭的角運動穩(wěn)定性，推導了彈箭的非線性角運動方程組，給出彈箭的非線性角運動Hopf分岔分析方法。以某型火箭彈高原試驗為例，選取空氣密度作為分岔參數(shù)，采用霍爾維茨判據(jù)判斷了系統(tǒng)的穩(wěn)定性，并確定了分岔點。由中心流形定理對系統(tǒng)進行降維，計算了Hopf分岔的3階規(guī)范形，并作出了系統(tǒng)的分岔圖，分析了分岔參數(shù)對極限環(huán)擺幅的影響。進行了仿真驗證，結果表明，采用分岔分析方法能準確判斷系統(tǒng)的穩(wěn)定性及分析系統(tǒng)的極限環(huán)運動。

兵器科學與技術；非線性角運動；運動穩(wěn)定性；Hopf分岔

0 引言

在火箭彈的高原遙測試驗過程中，觀測到大射角射擊時多次出現(xiàn)大攻角錐形運動，使得射程減小，這種現(xiàn)象可能是非線性運動造成的。對于彈箭的非線性角運動分析，文獻［1］中推導了復攻角模型，并借助平均法等定性方法理論分析了不同情況下的非線性角運動規(guī)律。文獻［2］在文獻［1］模型的基礎上研究了非旋轉大長徑比彈箭在非線性氣動力作用下產生極限平面擺動的機理及其抑制措施。文獻［3］推導一個較為復雜的保留了幾何非線性和氣動非線性的復攻角模型，但未涉及到具體分析。文獻［4］通過對大長徑比無控低旋火箭彈的卷弧尾翼在3種不同安裝位置時的氣動特性進行數(shù)值計算，得出火箭彈在飛行過程中確實會出現(xiàn)錐形運動且尾翼安裝方式對錐形運動會產生影響，證明了反裝反向滾轉卷弧形尾翼可以抑制低旋無控火箭彈的錐形運動。文獻［5］考察了章動運動條件下的自旋彈體運動模型，并分析了線性化模型的穩(wěn)定性，計算了自旋彈體的穩(wěn)定性臨界轉速與臨界錐角。文獻［6］建立了章動和進動的復合運動模型，根據(jù)線性化模型分析了彈體章動及進動運動的相互關系，給出了自旋彈體章動運動穩(wěn)定性條件。文獻［7］綜合考慮彈體姿態(tài)運動和位移運動建立了旋轉彈錐形運動的動力學模型，采用李雅普諾夫一級近似方法，給出了彈道頂點附近彈體錐形運動的穩(wěn)定判據(jù)。文獻［8］給出了火箭錐形運動的微分方程組，文獻［9-10］在該方程組的基礎上通過小偏差線性化方法研究了火箭彈在小錐角情況下圓錐運動漸進穩(wěn)定條件及收斂速度計算方法，并提出了提高火箭彈圓錐運動穩(wěn)定性的方法。文獻［11-12］基于陀螺線性擾動運動方程研究了旋轉導彈的錐形運動穩(wěn)定性問題。事實上，錐形運動對應于常微分方程中的極限環(huán)問題，因此，可采用常微分方程穩(wěn)定性分岔理論進行研究。

本文建立了新的彈箭非線性角運動方程組，該方程組盡可能完整地保留了幾何非線性及氣動力非線性。給出了基于中心流形定理和規(guī)范形理論的彈箭非線性角運動穩(wěn)定性Hopf分岔分析方法。結合某型火箭彈在高原射擊時的參數(shù)，分析了以密度作為分岔參數(shù)時火箭彈非線性角運動的分岔特性。最后在給定初始條件下進行了數(shù)值模擬，驗證了分岔分析方法結果正確性。

1 彈箭非線性角運動方程

1.1 坐標系的定義［1］

1）地面坐標系OExyz:其原點OE在炮口斷面中心，OEx軸沿著水平線指向目標，OEy軸垂直于OEx軸，鉛直向上為正，OEz軸按右手法則確定。地面坐標系為動坐標系，在忽略地球自轉的情況下可將其視為慣性系。

2）基準坐標系Oxryrzr:由地面坐標系平移至彈箭質心O而成，隨質心一起平動。

3）彈道坐標系Oxtytzt:其Oxt軸沿質心速度矢量v的方向，Oyt軸在包含速度矢量v的鉛垂面內垂直于Oxt軸，向上為正，Ozt軸按右手法則確定。

4）第2彈軸坐標系Oξηζ:其Oξ軸為彈軸，指向彈頭方向為正。該系由彈道坐標系經(jīng)過繞Ozt軸和Oη軸的分別旋轉δ1和δ2而來，則Oη軸在包含速度矢量v的鉛垂面內垂直于Oξ軸，向上為正。

該系還可由基準系經(jīng)三次旋轉而來，第一次是基準坐標系繞Ozn軸正向右旋φ1角到達Oξ′η′zn，第二次是Oξ′η′zn系繞Oη′軸負向右旋φ2角到達Oξη′ζ′，第三次是Oξη′ζ′系繞Oξ左旋β角得來。

第2彈軸坐標系的轉動角速度ω1為

5）彈體坐標系OXbYbZb:其OXb軸與彈軸平行；OYb軸垂直O(jiān)Xb軸，并在彈箭縱向平面內，向上為正；OZb軸滿足右手法則。

此坐標系可由第2彈軸坐標系繞彈軸轉過-β+γ角而得到，其轉動角速度矢量為

1.2 坐標轉換矩陣

1）彈道坐標系與第2彈軸坐標系間的轉換矩陣L（δ1，δ2）為

2）基準坐標系與彈道坐標系間的轉換矩陣L（θ1，ψ2）為

式中:θ1、ψ2分別為速度高低角和方位角。

1.3 彈箭非線性角運動方程推導

1.3.1 彈箭質心運動的動力學方程

彈箭質心運動方程組的矢量形式描述［1］如下:

式中:m為彈丸質量；F為彈丸受到的氣動力的合力；ωt為彈道坐標系的轉動角速度。將方程向第2彈軸坐標系投影得

式中:Fx、Fy、Fz為氣動力F在彈道坐標系3軸上的分量。（6）式等號左邊兩項計算如下:

將（7）式、（8）式代入（6）式，并在等號兩邊同時乘上LT（δ1，δ2），則有

式中:G為彈箭重力。

整理（9）式得

1.3.2 彈箭繞心運動的動力學方程

彈箭繞質心的轉動方程用矢量形式描述［1］如下:

式中:H為彈箭對質心的動量矩；M為彈箭所受外力對質心的合力矩。

方程左邊的兩項計算如下:

式中:A、C分別為彈箭的赤道轉動慣量和極轉動慣量；彈體坐標系的轉動角速度ω和第2彈軸坐標系的轉動角速度ω1的表達式分別為

整理（11）式～（14）式可以得到以第2彈軸坐標系為參考系的繞心運動的動力學方程組為

式中:Mξ、Mη、Mζ為力矩M在第2彈軸坐標系3軸上的分量。

1.3.3 彈箭非線性角運動方程

對于彈箭的角運動，主要關心的是δ1、δ2、ωη、ωζ的變化情況，故選取這4個量作為狀態(tài)變量x=［δ1δ2ωηωζ］T.

彈箭飛行過程中，有sin β≈sin ψ2sin δ1/cos δ2，δ2≈φ2-ψ2，則

角度單位為弧度時，ψ2和的量級為10-4和 10-5，在攻角和擺動角速度較大時，經(jīng)估算分析，相對于ωζtan φ2為小量，故取、φ2≈δ2.則ω和 ω1關系為

將第2彈軸坐標系的轉動角速度代入到δ1、δ2、ωη、ωζ的方程中，不計重力，得到彈箭的非線性角運動方程組如下:

根據(jù)文獻［1，13］，考慮氣動力和力矩的最簡非線性形式時，升力系數(shù)導數(shù)、靜力矩系數(shù)、赤道阻尼力矩系數(shù)、馬格努斯力矩系數(shù)可統(tǒng)一表示為以下多項式形式:

式中:δ為彈箭的總攻角，cos δ=cos δ2cos δ1；ct為上述的氣動系數(shù)；c0、c2分別為氣動系數(shù)的線性項和立方項。

方程組中包含了三角函數(shù)及其乘積，推導過程未涉及到三角函數(shù)的線性化問題，這體現(xiàn)了該方程組保留了幾何非線性。因此，該方程組具有普遍性，適用于彈箭的非線性角運動穩(wěn)定性分析。

2 彈箭非線性角運動分岔分析方法

2.1 非線性角運動方程分岔點計算

將彈箭所受的氣動力和力矩代入到角運動方程中得

式中:S、l、d分別為彈箭的參考面積、長度和直徑。

令方程組（20）式等號右邊為0，求得O［0000］T為系統(tǒng)的平衡點。（20）式在平衡點處的雅克比矩陣為

求出A的特征多項式，選取分岔參數(shù)μ，當參數(shù)μ變化時，根據(jù)霍爾維茨判據(jù)判斷系統(tǒng)穩(wěn)定性的變化，穩(wěn)定性發(fā)生改變的點即為分岔點μ0.此時特征多項式有一對純虛的特征根和一對實部為負的特征根（實部為正時，系統(tǒng)不穩(wěn)定，這是實際中不關心的），求出系統(tǒng)在分岔點處的特征值λ1、λ2、λ3、λ4及特征向量ξ1、ξ2、ξ3、ξ4.

2.2 非線性角運動方程降維

四維角運動方程直接進行分岔分析比較困難，需要對其進行降維。在分岔點處，系統(tǒng)特征多項式有一對純虛的特征根和一對實部為負的特征根，此時存在二維中心流形，可采用中心流形定理對系統(tǒng)進行降維［14］。

引入非奇異線性變換x=py，其中y=［y1y2y3y4］T，其中p為A的第2和第4個特征值所對應的特征向量的實部和虛部所構成的方陣。令u=［y1y2］T，v=［y3y4］T，代入系統(tǒng)方程得

式中:B、C分別為兩個方程的雅克比矩陣；F（u，v）和G（u，v）為方程組的非線性項。

由中心流形定理，可以把v表示為v=h（u），其中h（0）=h′（0）=0.為確定h（u），將v=h（u）代入（22）式的第2式，利用求導的鏈式法則，有

再利用第1式整理得到h（u）的微分方程

精確求解h（u）有困難，可以根據(jù)所需要中心流形的階數(shù)假設方程

將（25）式代入到h（u）的微分方程，根據(jù)各項系數(shù)對應相等解出aijk、bijk，即得到v=h（u）的表達式.將h（u）代入（22）式中的第1式，即可得到中心流形上流的約化方程

2.3 非線性角運動方程極限環(huán)計算

將約化后的方程寫成如下形式:

根據(jù)規(guī)范形理論［15］，其3階規(guī)范形為

引入極坐標變換y1=rcos θ，y2=rsin θ，得到如下方程:

由此得到系統(tǒng)的近似平衡方程cμr+ar3=0.系統(tǒng)的極限環(huán)由cμr+ar3=0的非零解決定。極限環(huán)的擺幅，周期T=2π/ω.

當a≠0，c≠0時，系統(tǒng)在μ=0處出現(xiàn)Hopf分岔，其分岔特性為

1）當c＞0和a＞0時，平衡點對μ＞0不穩(wěn)定，對μ＜0漸進穩(wěn)定，當μ＜0存在不穩(wěn)定極限環(huán)；

2）當c＞0和a＜0時，平衡點對μ＞0不穩(wěn)定，對μ＜0漸進穩(wěn)定，當μ＞0存在漸進穩(wěn)定極限環(huán)；

3）當c＜0和a＞0時，平衡點對μ＞0漸進穩(wěn)定，對μ＜0不穩(wěn)定，當μ＞0存在不穩(wěn)定極限環(huán)；

4）當c＜0和a＜0時，平衡點對μ＞0漸進穩(wěn)定，對μ＜0不穩(wěn)定，當μ＜0存在漸進穩(wěn)定極限環(huán)。

3 火箭彈非線性角運動計算分析

根據(jù)外彈道學的理論［13］，尾翼彈的飛行不穩(wěn)定主要是由非線性馬格努斯力矩造成的。高原大射角射擊時，火箭彈主動段后彈道上空氣密度明顯降低，導致馬格努斯力矩的影響加劇，從而出現(xiàn)飛行不穩(wěn)定。降低射角后射擊，飛行不穩(wěn)定現(xiàn)象明顯減少，這是因為彈道上空氣密度增加，有效抑制了馬格努斯效應的不利影響。因此，密度的變化對于馬格努斯效應的影響程度較大。本節(jié)主要分析考慮立方馬格努斯力矩后，選取空氣密度作為分岔參數(shù)時，某型火箭彈的非線性角運動穩(wěn)定性的情況。

3.1 角運動方程分岔點計算

將火箭彈參數(shù)和所受的力和力矩代入到非線性角運動方程中，求得系統(tǒng)在平衡點O處的雅克比矩陣

它的特征多項式可以寫成如下形式:

式中:a4=1；a3=3.0242ρ；a2=0.0875+1 324.944 2ρ+ 3.4063ρ2；a1=17.401 8ρ+2 003.468 8ρ2+1.693 3ρ3；a0=4.387 3×105ρ2+741.855 5ρ3+0.313 5ρ4.

根據(jù)四維系統(tǒng)的霍爾維茨判據(jù):穩(wěn)定條件是特征方程的所有系數(shù)為正數(shù)，還要Δ3＞0.可以看出角方程雅克比矩陣的特征方程所有系數(shù)都為正數(shù)，因此只需要判斷Δ3的符號。由

計算可得:當ρ＜0.5623時，Δ3＜0，系統(tǒng)平衡點O是不穩(wěn)定的；ρ＞0.562 3時，Δ3＞0，平衡點O是穩(wěn)定的。所以，為分岔點，此時高度約為7400 m.

3.2 角運動方程降維

將方程組中的三角函數(shù)在零點鄰域內進行泰勒展開，忽略三次以上的高階量，令，進行非奇異線性變換x=py，系統(tǒng)化成

根據(jù)中心流形定理，計算其2階中心流形

忽略3階以上的項，中心流形上的約化方程為

這樣就利用中心流形定理將原四維系統(tǒng)化成了二維系統(tǒng)。

3.3 角運動極限環(huán)計算

本節(jié)計算密度對極限環(huán)的影響，包括兩部分:極限環(huán)的產生條件及極限環(huán)幅值和周期的計算。

3.3.1 角運動極限環(huán)產生條件

根據(jù)規(guī)范形理論，約化后的系統(tǒng)3階規(guī)范形為

式中:c=-0.378 0；ω=19.454 0；e=-17.165 1；a=-0.002 9；周期T=2π/ω=0.32 s.

引入極坐標變換，得到系統(tǒng)的近似平衡方程0.378 0μr+0.002 9r3=0.其分岔圖如圖1所示，其中橫坐標μ為系統(tǒng)分岔參數(shù)，縱坐標r表示Hopf分岔。橫軸上實線代表穩(wěn)定的平衡點，虛線代表不穩(wěn)定平衡點。

圖1 系統(tǒng)Hopf分岔圖Fig.1 Hopf bifurcation diagram of system

由分岔圖可以看出:當μ＞0時，平衡方程的只有唯一解r=0，此時平衡點穩(wěn)定；當μ＜0時，平衡

3.3.2 角運動極限環(huán)計算

由二維中心流形

則y3、y4表示為

由x=py換回物理坐標得

空氣密度ρ分別取0.3 kg/m3、0.4 kg/m3、0.5 kg/m3時，系統(tǒng)的極限環(huán)如圖2所示。

由圖2中可以看出，密度對極限環(huán)的幅值影響比較大，當密度減小時，極限環(huán)幅值增大。當密度為0.3 kg/m3時，極限環(huán)幅值約為13°.

4 數(shù)值模擬

當空氣密度取ρ=1.2 kg/m3，初始條件為［-1.1° 1.1° 0.6 rad/s -0.6 rad/s］T時角運動如圖3所示。

圖2 不同密度時的角運動極限環(huán)Fig.2 Limit cycles of angular motion for different densities

圖3 ρ=1.2 kg/m3時角運動軌跡相圖Fig.3 Trajectory diagram of angular motion for ρ=1.2 kg/m3

當空氣密度取ρ=0.3 kg/m3，初始條件為［-8.0° 3.8° 7.5 rad/s -2.1 rad/s］T時角運動如圖4所示。

圖4 ρ=0.3 kg/m3時角運動軌跡相圖Fig.4 Trajectory diagram of angular motion for ρ=0.3 kg/m3

從圖3和圖4中可以看出:在初始值給定的情況下，ρ=1.2 kg/m3時，系統(tǒng)的運動為趨近于平衡點；ρ=0.3 kg/m3時，系統(tǒng)的運動為趨近于極限環(huán)；當ρ=0.3 kg/m3時，數(shù)值模擬得到δ1、δ2的極限環(huán)擺幅為13°，和計算相同。

5 結論

本文建立了彈箭非線性角運動方程，給出了彈箭非線性角運動Hopf分岔分析方法。選取某型火箭彈高原射擊時的數(shù)據(jù)，以密度作為分岔參數(shù)，計算了角運動的分岔特性。結果顯示，密度在0.5623 kg/m3時，角運動開始出現(xiàn)不穩(wěn)定，此時高度約為7 400 m，這與高原試驗中開始出現(xiàn)飛行不穩(wěn)定的高度基本相同。當密度小于0.562 3 kg/m3時，火箭彈在非線性馬格努斯力矩作用下角會出現(xiàn)極限環(huán)運動，且極限環(huán)的擺幅隨著密度的減小而增大。當密度為0.3 kg/m3時，極限環(huán)的擺幅約為13°.最后數(shù)值模擬驗證了分析方法的正確性。本文建立的非線性角運動模型及給出的非線性角運動Hopf分岔分析方法對彈箭的飛行理論、結構設計、試驗分析等具有一定的應用價值。

（

）

［1］ 韓子鵬.彈箭外彈道學［M］.北京:北京理工大學出版社，2008. HAN Zi-Peng.Exterior ballistics of shells and rockets［M］.Beijing:Beijing Institute of Technology Press，2008.（in Chinese）

［2］ 李臣明，劉怡昕.大長徑比遠程彈箭的極限平面擺動及其抑制［J］.系統(tǒng)與仿真學報，2009，21（23）:7390-7393. LI Chen-ming，LIU Yi-xin.Limit plane swing motion and its restraining measure of un-rotary long-range missile with large ratio of length to diameter［J］.Journal of System Simulation，2009，21（23）:7390-7393.（in Chinese）

［3］ 徐明友.高等外彈道學［M］.北京:高等教育出版社，2003. XU Ming-you.Advanced exterior ballistics［M］.Beijing:Higher Education Press，2003.（in Chinese）

［4］ 王華畢，吳甲生.火箭彈錐形運動的數(shù)學仿真與抑制措施［J］.北京理工大學學報，2007，27（3）:196-199. WANG Hua-bi，WU Jia-sheng.Coning motion of rockets，its numerical simulation and restraint［J］.Transactions of Beijing Institute of Technology，2007，27（3）:196-199.（in Chinese）

［5］ 王華畢，吳甲生.火箭彈錐形運動穩(wěn)定性分析［J］.兵工學報，2008，29（5）:562-566. WANG Hua-bi，WU Jia-sheng.The coning motion stability analysis of rocket［J］.Acta Armamentarii，2008，29（5）:562-566.（in Chinese）

［6］ 閆曉勇，楊樹興，張成.基于章動運動理論的火箭彈錐形運動穩(wěn)定性分析［J］.兵工學報，2009，30（10）:1291-1296. YAN Xiao-yong，YANG Shu-xing，ZHANG Cheng.Analysis of stability for coning motion of rockets based on theory of nutation movement［J］.Acta Armamentarii，2009，30（10）:1291-1296.（in Chinese）

［7］ 李克勇，趙良玉，周偉.一類旋轉彈在高空中的錐形運動穩(wěn)定性［J］.動力學與控制學報，2012，10（3）:239-243. LI Ke-yong，ZHAO Liang-yu，ZHOU Wei.Stability for coning motion of a spinning projectile［J］.Journal of Dynamic and Control，2012，10（3）:239～243.（in Chinese）

［8］ Mao X R，Yang S X，Xu Y.Research on the coning motion of wrap around fin projectiles［J］.Canadian Aeronautics and Space Journal，2006，52（3）:119-125.

［9］ 趙良玉，楊樹興.卷弧翼火箭彈圓錐運動收斂速度計算方法［J］.固體火箭技術，2009，32（1）:15-19. ZHAO Liang-yu，YANG Shu-xing.Research on convergence speed of coning motion of wrap around fin rockets［J］.Journal of Solid Rocket Technology，2009，32（1）:15-19.（in Chinese）

［10］ 趙良玉，楊樹興，焦清介.提高卷弧翼火箭彈圓錐運動漸進穩(wěn)定性的幾個方法［J］.固體火箭技術，2010，33（4）:369-372. ZHAO Liang-yu，YANG Shu-xing，JIAO Qing-jie.Several methods for improving asymptotic stability of coning motion of wrap around fin rockets［J］.Journal of Solid Rocket Technology，2010，33（4）:369-372.（in Chinese）

［11］ 任天榮，馬建敏.基于陀螺力學的旋轉導彈錐形運動分析［J］.宇航學報，2010，31（9）:2082-2087. REN Tian-rong，MA Jian-min.Coning motion analysis of spinning missile based on gyro dynamics［J］.Journal of Astronautics，2010，31（9）:2082-2087.（in Chinese）

［12］ 任天榮，馬建敏.旋轉彈錐形運動發(fā)生區(qū)間及頻率特性研究［J］.固體火箭技術，2014，37（3）:295-300. REN Tian-rong，MA Jian-min.Research on activating region and frequency characteristics of coning motion for spinning missiles［J］.Journal of Solid Rocket Technology，2014，37（3）:295-300.（in Chinese）

［13］ McCoy R L.Modern exterior ballistics［M］.Atglen，Pennsylvania:Schiffer Publishing Ltd，1999.

［14］ 張琪昌，王洪禮，竺致文，等.分岔與混沌理論及應用［M］.天津:天津大學出版社，2005. ZHANG Qi-chang，WANG Hong-li，ZHU Zhi-wen，et al.Bifurcation and chaos theory and its application［M］.Tianjin:Tianjin University Press，2005.（in Chinese）

［15］ 陸啟韶.常微分方程與動力系統(tǒng)［M］.北京:北京航空航天大學出版社，2010. LU Qi-shao.Ordinary differential equations and dynamical systems［M］.Beijing:Beijing University of Aeronautics and Astronautics Press，2010.（in Chinese）

Hopf Bifurcation Analysis of Nonlinear Angular Motion Stability of Projectile

ZHONG Yang-wei，WANG Liang-ming，F(xiàn)U Jian，CHANG Si-jiang
（School of Energy and Power Engineering，Nanjing University of Science and Technology，Nanjing 210094，Jiangsu，China）

In order to analyze the angular motion stability of projectile，the equations of the nonlinear angular motion are derived，and the Hopf bifurcation analysis method of the nonlinear angular motion of projectile is given.Taking a rocket plateau test as an example，the air density is selected as the bifurcation parameter，and the Hurwitz criterion is used to judge the stability of the system.The bifurcation point is determined.Center manifold theory is proposed to reduce the system dimension，and then a three-order normal form of Hopf bifurcation is calculated by plotting the bifurcation diagram.In addition，the effect of the bifurcation parameter on the swing of the limit cycle is analyzed.The numerical simulations show that the bifurcation analysis method can be used to judge the stability of the system correctly and analyze the motion of limit cycle in the system accurately.

ordnance science and technology；nonlinear angular motion；motion stability；Hopf bifurcation

TJ714

A

1000-1093（2015）07-1195-08

10.3969/j.issn.1000-1093.2015.07.007

2014-05-08

國家自然科學基金項目（11402117）

鐘揚威（1989—），男，博士研究生。E-mail:zyw_601@163.com；王良明（1963—），男，教授，博士生導師。E-mail:lmwang802@163.com