彈箭非線性角運(yùn)動(dòng)周期解穩(wěn)定性分析

鐘揚(yáng)威，王良明，常思江，傅 健

（南京理工大學(xué) 能源與動(dòng)力工程學(xué)院，南京210094）

在火箭彈的高原遙測(cè)試驗(yàn)中，觀測(cè)到在大射角射擊時(shí)出現(xiàn)多次近彈現(xiàn)象，姿態(tài)測(cè)量結(jié)果表明火箭彈在主動(dòng)段結(jié)束時(shí)出現(xiàn)約20°左右的大攻角錐形運(yùn)動(dòng)。查閱資料發(fā)現(xiàn)，美國(guó)在進(jìn)行Nitehawk探空火箭飛行試驗(yàn)時(shí)，50多次飛行試驗(yàn)出現(xiàn)了20多次發(fā)散的錐形運(yùn)動(dòng)。美國(guó)的69.85mm機(jī)載火箭彈在亞音速及超音速風(fēng)洞實(shí)驗(yàn)時(shí)都曾出現(xiàn)發(fā)散的錐形運(yùn)動(dòng)［1］。西班牙的Ternel 140mm火箭彈28次飛行試驗(yàn)中出現(xiàn)了9次錐形運(yùn)動(dòng)，使飛行速度在1.5s內(nèi)降低了60%［2］。尾翼穩(wěn)定火箭彈的這些角運(yùn)動(dòng)規(guī)律用線性理論難以解釋，可能是非線性運(yùn)動(dòng)造成的。

目前，許多學(xué)者對(duì)尾翼穩(wěn)定旋轉(zhuǎn)彈箭的錐形運(yùn)動(dòng)穩(wěn)定性進(jìn)行了研究。文獻(xiàn)［3］研究了卷弧尾翼安裝位置對(duì)大長(zhǎng)徑比無(wú)控低旋火箭彈錐形運(yùn)動(dòng)的影響，證明了通過(guò)反裝反向卷弧形尾翼可以抑制火箭彈的錐形運(yùn)動(dòng)。文獻(xiàn)［4］建立了無(wú)控低旋火箭彈的錐形運(yùn)動(dòng)方程組，并分析了側(cè)向力矩對(duì)錐形運(yùn)動(dòng)穩(wěn)定性的影響。文獻(xiàn)［5］建立了彈體錐形運(yùn)動(dòng)的模型，研究了章動(dòng)運(yùn)動(dòng)對(duì)錐形運(yùn)動(dòng)極限環(huán)的影響。文獻(xiàn)［6］研究了大攻角情況下非對(duì)稱赤道阻尼力矩對(duì)遠(yuǎn)程火箭極限圓錐運(yùn)動(dòng)的影響。文獻(xiàn)［7－8］基于陀螺線性擾動(dòng)運(yùn)動(dòng)方程研究了旋轉(zhuǎn)導(dǎo)彈的錐形運(yùn)動(dòng)穩(wěn)定性問(wèn)題。文獻(xiàn)［9］對(duì)一類低速旋轉(zhuǎn)彈在高空且考慮非線性氣動(dòng)力時(shí)的錐形運(yùn)動(dòng)穩(wěn)定性進(jìn)行了分析。

事實(shí)上，彈箭的錐形運(yùn)動(dòng)對(duì)應(yīng)于非線性動(dòng)力學(xué)中的周期運(yùn)動(dòng)，本文采用非線性動(dòng)力學(xué)中的周期解分岔及穩(wěn)定性分析等相關(guān)理論，基于新的彈箭非線性角運(yùn)動(dòng)方程，提出了分析和計(jì)算彈箭非線性角運(yùn)動(dòng)周期解的方法，為分析彈箭非線性角運(yùn)動(dòng)規(guī)律提供了理論基礎(chǔ)。

1 彈箭非線性角運(yùn)動(dòng)方程組

對(duì)于彈箭的角運(yùn)動(dòng)，主要關(guān)心的是攻角δ1，δ2及擺動(dòng)角速度ωη，ωζ的變化情況。故選取這4個(gè)量作為狀態(tài)變量x＝（δ1δ2ωηωζ）T，根據(jù)文獻(xiàn)［10］中的彈箭運(yùn)動(dòng)方程組，通過(guò)推導(dǎo)得到彈箭的非線性角運(yùn)動(dòng)方程組：

式中：Fy，F(xiàn)z分別為彈箭所受的總氣動(dòng)力在速度坐標(biāo)系相應(yīng)軸上的投影；m為彈箭的質(zhì)量；v為彈箭的速度；Mη和Mζ為彈箭所受的合力矩在彈軸坐標(biāo)系相應(yīng)軸上的投影；ωξ為彈箭的滾轉(zhuǎn)角速度；JA，JC分別為彈箭的赤道轉(zhuǎn)動(dòng)慣量和極轉(zhuǎn)動(dòng)慣量。氣動(dòng)力的計(jì)算公式見(jiàn)文獻(xiàn)［11］。

將角運(yùn)動(dòng)方程寫(xiě)成規(guī)范形式，可以表示為＝f（x，μ），其中x∈R4是狀態(tài)向量；μ∈Rm是分岔參數(shù)，是可以表示空氣密度、馬赫數(shù)、彈箭自轉(zhuǎn)速度、氣動(dòng)力和力矩系數(shù)等影響彈箭角運(yùn)動(dòng)的參數(shù)。當(dāng)選取的參數(shù)μ連續(xù)變化時(shí)，系統(tǒng)的拓?fù)浣Y(jié)構(gòu)在μ＝μ0處發(fā)生突然變化，則稱系統(tǒng)在μ＝μ0處出現(xiàn)分岔，并稱μ0為一個(gè)分岔值。彈箭非線性角運(yùn)動(dòng)周期解穩(wěn)定性分析主要解決兩方面的問(wèn)題：①?gòu)椉沁\(yùn)動(dòng)是否存在以及何時(shí)存在周期解；②周期解的幅值及周期計(jì)算和穩(wěn)定性分析。

2 角運(yùn)動(dòng)分岔序列計(jì)算

彈箭角運(yùn)動(dòng)方程是一個(gè)連續(xù)性方程，可以通過(guò)計(jì)算其Poincare映射將其離散化來(lái)分析其分岔行為［12］。Poincare映射是研究周期運(yùn)動(dòng)及其分岔的幾何方法，它將非線性系統(tǒng)轉(zhuǎn)化為與之本質(zhì)上完全等價(jià)而研究起來(lái)更為簡(jiǎn)單方便的點(diǎn)映射系統(tǒng)。Poincare映射的定義如下。

設(shè)Γ為系統(tǒng)的周期解，其最小周期為T，過(guò)Γ上任意一點(diǎn)x＊作與Γ橫截相交的n－1維超曲面。根據(jù)流形φt的連續(xù)性，存在x＊的某個(gè)鄰域U?Σ，使得U內(nèi)任何一點(diǎn)出發(fā)的軌道都可以在約T的時(shí)間內(nèi)再次回到Σ上，這樣就定義了x＊的一個(gè)鄰域U?Σ到Σ的映射P，P稱為系統(tǒng)的Poincare映射。

通過(guò)Poincare映射可以計(jì)算彈箭非線性角運(yùn)動(dòng)的分岔序列，其數(shù)值算法描述如下。

①選擇超曲面Σ。

彈箭角運(yùn)動(dòng)方程為非線性自治方程，自治系統(tǒng)沒(méi)有特定的方法選取Poincare截面，需要根據(jù)具體的運(yùn)動(dòng)情況進(jìn)行選取。本文選取彈箭的擺動(dòng)角速度ωη的導(dǎo)數(shù)通過(guò)0時(shí)的點(diǎn)作為軌道與Poincare截面的交點(diǎn)，即定義

②計(jì)算軌道與超曲面的交點(diǎn)。

選取分岔參數(shù)μ為某一確定值，對(duì)角運(yùn)動(dòng)方程進(jìn)行積分，計(jì)算經(jīng)過(guò)每一個(gè)步長(zhǎng)時(shí)的。重復(fù)以上步驟，直到連續(xù)的2個(gè)點(diǎn)｜t＝t1和｜t＝t2位于Σ的兩側(cè)，即有相反的符號(hào)，則可確定交點(diǎn)位于t1和t2之間的某一點(diǎn)。不斷平分區(qū)間［t1，t2］直至達(dá)到要求的精度。重復(fù)上述過(guò)程，可求出所給定的任意時(shí)間長(zhǎng)度內(nèi)的所有交點(diǎn)。

③計(jì)算分岔序列。

連續(xù)地改變分岔參數(shù)μ，重復(fù)步驟②，直到求出所要求的分岔區(qū)間的Poincare映射圖。將分岔區(qū)間內(nèi)的Poincare映射圖疊加即可得到分岔序列。

3 周期解計(jì)算及穩(wěn)定性分析

通過(guò)計(jì)算角運(yùn)動(dòng)的分岔序列得出系統(tǒng)存在周期解，如何迅速求出其周期解的幅值和周期是一個(gè)問(wèn)題。采用龍格庫(kù)塔法等數(shù)值積分方法通過(guò)積分來(lái)求周期解具有一定的盲目性且只能獲得穩(wěn)定的周期解。因此，在周期解的數(shù)值計(jì)算時(shí)一般不采用直接積分的方法。打靶法是一種常用的求解系統(tǒng)周期解的方法，本文采用推廣的打靶法［13－14］獲得彈箭非線性角運(yùn)動(dòng)周期解的幅值及周期，并通過(guò)Floquet理論分析周期解的穩(wěn)定性。

3.1 角運(yùn)動(dòng)周期解數(shù)值計(jì)算

彈箭角運(yùn)動(dòng)的周期是未知的，因此需要將周期作為一個(gè)參數(shù)一起參與打靶，從而求得角運(yùn)動(dòng)的周期解及周期?；谕茝V的打靶法求取彈箭非線性角運(yùn)動(dòng)周期解的步驟如下。

①引入變量τ，變換時(shí)間尺度t＝Tτ，角運(yùn)動(dòng)方程變?yōu)?/p>

②選取初始迭代向量x0＝從τ＝0到τ＝1對(duì)式（3）進(jìn)行積分，得到系統(tǒng)在周期末時(shí)刻的值x1＝），并計(jì)算誤差向量r，其中：

式中：i，j＝1，…，4；xi，xj分別表示為Kronecker符號(hào)。

⑤將步驟④中所求的值代入式（5）中，求得?ri／，?ri／?T0的值。

⑥將ri＝0在xi和T0鄰域內(nèi)作泰勒展開(kāi)，忽略高階項(xiàng)，生成的線性方程組為

⑧檢查步驟⑦求出的增量值是否滿足所需的精度要求，若滿足則結(jié)束迭代，得到角運(yùn)動(dòng)的周期解和周期；否則令＝＋Δxi和T0＝T0＋ΔT轉(zhuǎn)第②步繼續(xù)進(jìn)行迭代計(jì)算。

3.2 周期解穩(wěn)定性分析方法

Floquet理論是從攝動(dòng)方程零解穩(wěn)定性判別未擾動(dòng)周期運(yùn)動(dòng)穩(wěn)定性的思想得到的穩(wěn)定性理論，利用它可以分析非線性角運(yùn)動(dòng)周期解的穩(wěn)定性問(wèn)題［15］。

對(duì)于非線性系統(tǒng)的攝動(dòng)方程：

式中：A（t＋T）＝A（t）是一個(gè)周期為T的矩陣函數(shù)。若x（t）為方程（7）的一個(gè)基本解矩陣，根據(jù)Floquet理論，則必存在一個(gè)非奇異的T周期矩陣φ（t）＝φ（t＋T）和一個(gè)常矩陣D，使得x（t）＝φ（t）etD。因?yàn)椋╰）＝A（t）x（t），則有：

從而

式中：eTD＝C為一常數(shù)矩陣，稱為系統(tǒng)的離散狀態(tài)傳遞矩陣，C的特征值λ稱為Floquet乘子。根據(jù)定義，狀態(tài)轉(zhuǎn)移矩陣C可以取為3.1中步驟⑦所計(jì)算出的系數(shù)矩陣。

可以利用主導(dǎo)Floquet乘子（絕對(duì)值最大的Floquet乘子）得到系統(tǒng)周期解的穩(wěn)定性判據(jù)：

①當(dāng)λ的絕對(duì)值小于1時(shí)，周期解穩(wěn)定；

②當(dāng)λ的絕對(duì)值大于1時(shí)，周期解不穩(wěn)定；

③當(dāng)λ的絕對(duì)值等于1時(shí)，周期解臨界穩(wěn)定。

4 數(shù)值仿真與分析

高原射擊時(shí)，火箭彈主動(dòng)段結(jié)束后其所處的飛行空域較高，空氣密度較低，在持續(xù)重力及較小的氣動(dòng)力和氣動(dòng)力矩作用下，容易誘導(dǎo)產(chǎn)生較大的攻角［9］。大攻角飛行時(shí)，彈箭的氣動(dòng)系數(shù)存在明顯的非線性，需要考慮氣動(dòng)系數(shù)的非線性項(xiàng)。

文獻(xiàn)［11］通過(guò)閃光攝影和飛行試驗(yàn)獲得了一些彈箭的非線性馬格努斯力矩系數(shù)，給出了非線性馬格努斯力矩的流體動(dòng)力學(xué)機(jī)理解釋，并指出了彈箭的非線性馬格努斯力矩與極限環(huán)的運(yùn)動(dòng)關(guān)系。筆者所在的課題組在進(jìn)行氣動(dòng)力數(shù)值計(jì)算時(shí)，發(fā)現(xiàn)在較大攻角時(shí)火箭彈的馬格努斯力矩系數(shù)存在嚴(yán)重非線性，它隨馬赫數(shù)和攻角的變化較為明顯。因此，高原彈箭射擊試驗(yàn)中出現(xiàn)的大攻角錐形運(yùn)動(dòng)可能是在高空低密度情況下由非線性馬格努斯力矩導(dǎo)致的。

根據(jù)文獻(xiàn)［11］，考慮線性馬格努斯力矩系數(shù)CMpα0和立方馬格努斯力矩系數(shù)CMpα2時(shí)，馬格努斯力矩系數(shù)CMpα可表示為以下的多項(xiàng)式形式：

式中：sin2αt＝（sinδ1cosδ2）2＋sin2δ2，αt為總攻角。

以某型火箭彈高原射擊時(shí)一些參數(shù)為主要仿真條件：①?gòu)椡栀|(zhì)量45kg，極轉(zhuǎn)動(dòng)慣量0.115kg·m2，赤道轉(zhuǎn)動(dòng)慣量30.63kg·m2，彈丸直徑0.122m，彈丸長(zhǎng)度2.8m；②陣地海拔為4 000m，火箭彈主動(dòng)段結(jié)束時(shí)刻飛行速度為1 012.3m／s，彈丸軸向角速度為79.006 4rad／s；③線性升力系數(shù)為6.5，線性靜力矩系數(shù)為－1.254 6，赤道阻尼力矩系數(shù)為0.46；④空氣密度取0.562 3kg／m3。

取CMpα0＝8，立方馬格努斯力矩系數(shù)CMpα2分別取10，20，30，40。仿真得到的角運(yùn)動(dòng)相圖和Poincare截面圖如圖1所示，圖中的＊表示龐加萊截面上的點(diǎn)。

圖1 角運(yùn)動(dòng)相圖

從圖1可以看出，CMpα2為正時(shí)，龐加萊截面上存在2個(gè)點(diǎn)，角運(yùn)動(dòng)存在周期解。進(jìn)一步計(jì)算，CMpα2為負(fù)時(shí)，龐加萊截面上沒(méi)有點(diǎn)，角運(yùn)動(dòng)發(fā)散。

取CMpα2＝30，取線性馬格努斯力矩系數(shù)CMpα0由－10～0變化，仿真得到的Poincare映射分岔圖，如圖2所示。

圖2 角運(yùn)動(dòng)分岔圖

由分岔圖看出，當(dāng)CMpα0＞－5.8時(shí)，角運(yùn)動(dòng)漸進(jìn)穩(wěn)定于唯一的零平衡位置；當(dāng)CMpα0＜－5.8時(shí)，零平衡位置不穩(wěn)定，分岔出周期解，角運(yùn)動(dòng)發(fā)生Hopf分 岔，CMpα0＝ －5.8 稱 為 分 岔 點(diǎn)。CMpα0＜－5.8后，角運(yùn)動(dòng)處于周期運(yùn)動(dòng)，周期運(yùn)動(dòng)的幅值隨CMpα0的減小而增大。線性馬格努斯力矩系數(shù)CMpα0分別?。?，－7，－8，－9時(shí)，通過(guò)推廣的打靶法計(jì)算得到的攻角運(yùn)動(dòng)周期解如圖3所示。由圖3可以看出，線性馬格努斯力矩系數(shù)CMpα0對(duì)角運(yùn)動(dòng)周期解的幅值影響比較大，CMpα0減小時(shí)，周期解幅值增大，與圖2得出的結(jié)論一致。當(dāng)線性馬格努斯力矩為－9時(shí)，攻角運(yùn)動(dòng)周期解的幅值約為19°。

圖3 角運(yùn)動(dòng)周期解

攻角運(yùn)動(dòng)周期解的周期T及主導(dǎo)Floquet乘子λ如表1所示。

表1 周期T及主導(dǎo)Floquet乘子λ

由表1可以看出，角運(yùn)動(dòng)周期解的周期隨線性馬格努斯力矩系數(shù)CMpα0的變化較小，且主導(dǎo)Floquet乘子λ的絕對(duì)值都小于1，說(shuō)明周期解是穩(wěn)定的。

5 結(jié)束語(yǔ)

本文基于一般形式的彈箭非線性角運(yùn)動(dòng)微分方程組，提出了基于Poincare映射的彈箭非線性角運(yùn)動(dòng)分岔序列的數(shù)值計(jì)算方法，給出了通過(guò)推廣的打靶法計(jì)算彈箭角運(yùn)動(dòng)周期解的幅值和周期的計(jì)算方法，并結(jié)合Floquet理論給出了周期解穩(wěn)定性的分析方法。以某型火箭彈為例，計(jì)算分析了立方馬格努斯系數(shù)作用下的角運(yùn)動(dòng)相圖及以線性馬格努斯力矩系數(shù)作為分岔參數(shù)時(shí)的非線性角運(yùn)動(dòng)規(guī)律。結(jié)果表明，在高空低密度的情況下，當(dāng)馬格努斯力矩系數(shù)達(dá)到一定的范圍時(shí)，火箭彈會(huì)出現(xiàn)穩(wěn)定的周期運(yùn)動(dòng)。

［1］李臣明，劉怡昕.轉(zhuǎn)速閉鎖對(duì)遠(yuǎn)程彈箭的影響［J］.彈道學(xué)報(bào)，2010，22（1）：45－48.LI Chen－ming，LIU Yi－xin.Effect of roll locking on long－range missile［J］.Journal of Ballistics，2010，22（1）：45－48.（in Chinese）

［2］MOROTE J，LIANO G.Stability analysis and flight trials of a clipped wrap around fin configuration［C］／／AIAA Atmospheric Flight Mechanics Conference and Exhibit.Providence，Rhode Island：AIAA，2004.

［3］王華畢，吳甲生.火箭彈錐形運(yùn)動(dòng)的數(shù)學(xué)仿真與抑制措施［J］.北京理工大學(xué)學(xué)報(bào)，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）

［4］王華畢，吳甲生.火箭彈錐形運(yùn)動(dòng)穩(wěn)定性分析［J］.兵工學(xué)報(bào)，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）

［5］閆曉勇，楊樹(shù)興，張成.基于章動(dòng)運(yùn)動(dòng)理論的火箭彈錐形運(yùn)動(dòng)穩(wěn)定性分析［J］.兵工學(xué)報(bào)，2009，30（10）：1 291－12 96.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）：1 291－1 296.（in Chinese）

［6］李臣明，劉怡昕.非對(duì)稱赤道阻尼力矩對(duì)遠(yuǎn)程火箭彈道的影響［J］.彈道學(xué)報(bào)，2009，21（2）：36－39.LI Chen－ming，LIU Yi－xin.Effect of roll locking on long－range missile［J］.Journal of Ballistics，2009，21（2）：36－39.（in Chinese）

［7］任天榮，馬建敏.基于陀螺力學(xué)的旋轉(zhuǎn)導(dǎo)彈錐形運(yùn)動(dòng)分析［J］.宇航學(xué)報(bào)，2010，31（9）：2 082－2 087.REN Tian－rong，MA Jian－min.Coning motion analysis of spinning missile based on gyro dynamics［J］.Journal of Astronautics，2010，31（9）：2 082－2 087.（in Chinese）

［8］任天榮，馬建敏.旋轉(zhuǎn)彈錐形運(yùn)動(dòng)發(fā)生區(qū)間及頻率特性研究［J］.固體火箭技術(shù)，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）

［9］李克勇，趙良玉，周偉.一類旋轉(zhuǎn)彈在高空中的錐形運(yùn)動(dòng)穩(wěn)定性［J］.動(dòng)力學(xué)與控制學(xué)報(bào)，2012，10（3）：239－243.LI Ke－yong，ZHAO Liang－yu，ZHOU Wei.Stability of conig motion of spinning rocket projectiles at high altitude［J］.Journal of Dynamic and Control，2012，10（3）：239－243.（in Chinese）

［10］韓子鵬.彈箭外彈道學(xué)［M］.北京：北京理工大學(xué)出版社，2014.HAN Zi－peng.Exterior ballistics of projectiles and rockets［M］.Beijing：Beijing Institute of Technology Press，2014.（in Chinese）

［11］ROBERT L.Modern exterior ballistics［M］.Atglen，Pennsylvania：Schiffer Publishing Ltd，1999：273－298.

［12］吳志強(qiáng)，張建偉.二元機(jī)翼極限環(huán)顫振復(fù)雜分岔［J］.工程力學(xué)，2008，25（2）：52－56.WU Zhi－qiang，ZHANG Jian－wei.Complicated bifurcation in limit cycle flutter of two dimensional airfoil［J］.Engineering Mechanics，2008，25（2）：52－56.（in Chinese）

［13］ROOSE D，LUST K，CHAMPENEYS A.A Newton－Picard shooting method for computing periodic solutions of large－scale dynamical systems［J］.Chaos，Solitons ＆ Fractals，1995，5（10）：1 913－1 925.

［14］李德信，徐健學(xué).求解非線性動(dòng)力系統(tǒng)周期解推廣的打靶法［J］.應(yīng)用力學(xué)學(xué)報(bào)，2003，20（4）：80－86.LI De－xin，XU Jian－xue.The generalized shooting method for solving the periodic solutions of nonlinear dynamical system［J］.Journal of Applied Mechanics，2003，20（4）：80－86.（in Chinese）

［15］張家忠.非線性動(dòng)力系統(tǒng)的運(yùn)動(dòng)穩(wěn)定性、分岔理論及其應(yīng)用［M］.西安：西安交通大學(xué)出版社，2010.ZHANG Jia－zhong.Motion stability，bifurcation theory and its application of nonlinear dynamical systems［M］.Xi’an：Xi’an Jiao Tong University Press，2010.（in Chinese）