楊志偉，王良明，鐘揚威,2，王垚,，張喜峰

(1.南京理工大學 能源與動力工程學院, 江蘇 南京 210094; 2.中國航天科工集團有限公司 第九總體設計部, 湖北 武漢 430040;3.北方華安工業(yè)集團有限公司, 黑龍江 齊齊哈爾 161046)

0 引言

彈丸的穩(wěn)定性一直是外彈道研究的重點和難點。早在20世紀60年代，美國外彈道專家Murphy就采用復攻角方法建立了彈丸的線性角運動方程，很好地描述了彈丸的角運動特性[1]。該方法隨后也被眾多國內(nèi)外外彈道學者[2-5]所采納，解決了早期彈丸穩(wěn)定性方面的問題。

隨著武器工業(yè)的發(fā)展，彈丸的氣動外形和使用環(huán)境都在不斷地發(fā)生變化，從而出現(xiàn)了一些新的彈丸角運動現(xiàn)象，如西班牙140 mm火箭彈實驗過程中的“近彈”現(xiàn)象、高原火箭彈飛行不穩(wěn)定現(xiàn)象。針對這些問題，國內(nèi)外學者進行了大量研究。文獻[4]從復攻角方程出發(fā)，采用振幅平面法分析了彈丸的非線性角運動特性，給出了非線性運動的動態(tài)穩(wěn)定性判據(jù)。文獻[6]在文獻[4]的基礎上，考慮幅值變化對振幅平面方程中根號項取值的影響，推導出了旋轉(zhuǎn)尾翼彈箭穩(wěn)定極限圓錐運動的解析判據(jù)。雖然這些方法能夠給出解析形式的穩(wěn)定性判據(jù)，但是在方程建立和求解過程中應用了大量假設和近似，削弱了方程的非線性。馬國梁等[7]通過線性化方法推導了彈丸角運動的狀態(tài)空間方程，提出了一種新的彈丸動態(tài)穩(wěn)定因子，并利用該穩(wěn)定因子解釋了高原環(huán)境對彈丸動態(tài)穩(wěn)定性的影響。鐘揚威等[8]針對火箭彈在高原條件下出現(xiàn)“近彈”問題，以密度作為分岔參數(shù)，利用中心流形定理在分岔點對彈丸角運動方程進行降維，然后對火箭彈的角運動進行Hopf分岔分析。文獻[9]采用相同的方法對一種內(nèi)置質(zhì)量塊的尾翼式彈道修正彈進行了Hopf分岔特性分析。文獻[8-9]采用數(shù)值方法求解彈箭非線性角運動方程，雖然無法得到解析解，但方程建立和求解過程假設很少，很大程度地保留了運動的非線性，針對具體問題所獲得解的精度非常高。

上述研究大多數(shù)是針對尾翼穩(wěn)定彈箭的非線性問題，對于旋轉(zhuǎn)穩(wěn)定彈丸很少涉及。相較于尾翼穩(wěn)定彈箭，旋轉(zhuǎn)穩(wěn)定彈丸需要通過高速自轉(zhuǎn)來維持其穩(wěn)定性，使得旋轉(zhuǎn)穩(wěn)定彈丸的動力學更加復雜。一方面，高速旋轉(zhuǎn)會產(chǎn)生陀螺效應，使得彈箭在俯仰、偏航和滾轉(zhuǎn)3個方向上的角運動高度耦合。另一方面，高速旋轉(zhuǎn)會產(chǎn)生非線性很強的馬格努斯力和力矩，增大了彈箭運動非線性分析的難度。文獻[10]提出，旋轉(zhuǎn)穩(wěn)定彈丸出現(xiàn)極限圓運動的前提條件是彈丸不滿足陀螺穩(wěn)定性。然而在近年來的射擊實驗中發(fā)現(xiàn)，對于已經(jīng)滿足了陀螺穩(wěn)定性和動態(tài)穩(wěn)定性條件的旋轉(zhuǎn)穩(wěn)定彈丸也會出現(xiàn)某些射擊條件下的運動失穩(wěn)。例如一些移動發(fā)射平臺(如艦炮、坦克等)，在其前進方向一側(cè)發(fā)射的彈丸運動穩(wěn)定，而另一側(cè)則不穩(wěn)定。又如美國M549彈丸在寒冷的冬季、2號裝藥發(fā)射條件下總是出現(xiàn)運動不穩(wěn)定現(xiàn)象。近年來頻發(fā)的旋轉(zhuǎn)彈丸“彈道炸”事故也集中發(fā)生在某些特定的射角和裝藥號下。

上述現(xiàn)象都與彈丸的初始條件相關，而初始條件對系統(tǒng)非線性穩(wěn)定性的影響可以轉(zhuǎn)化為計算平衡點的吸引域來分析。文獻[11]指出，系統(tǒng)平衡點的吸引域可用作非線性系統(tǒng)穩(wěn)定性相對危險性的度量。精確計算非線性系統(tǒng)的吸引域是一個非常復雜的問題，文獻[12-13]使用多項式平方和優(yōu)化方法進行計算。該方法只適用于多項式矢量所描述的動力學系統(tǒng)，并不適用于彈丸的角運動系統(tǒng)。

本文從旋轉(zhuǎn)彈丸的非線性角運動方程出發(fā)，給出一類平衡點吸引域邊界的計算方法。通過研究不同參數(shù)對此類平衡點吸引域邊界的影響，分析初始條件對彈丸穩(wěn)定性的影響，研究結(jié)果可為旋轉(zhuǎn)彈丸的氣動設計、射擊選擇提供參考。

1 旋轉(zhuǎn)穩(wěn)定彈丸動力學建模

1.1 坐標系和角度的定義

1)基準坐標系Oxnynzn：坐標原點O為彈丸質(zhì)心，Oxn軸沿水平線指向射擊方向，Oyn軸鉛直向上，Ozn由右手定則確定。

2)彈道坐標系Oxbybzb：坐標原點O為彈丸質(zhì)心，Oxb軸沿速度矢量方向，Oyb軸垂直于速度方向，Ozb由右手定則確定。速度矢量相對于基準坐標系的方位角可由速度高低角θa和速度方向角θd確定。

3)第一彈軸系Oξηfζf：坐標原點O為彈丸質(zhì)心，Oξ軸為彈軸，Oηf軸垂直于彈軸，Oζf由右手定則確定。彈軸相對于基準坐標系的方位角可由彈軸高低角φa和彈軸方向角φd確定。

1.2 旋轉(zhuǎn)穩(wěn)定彈丸非線性角運動模型

彈道坐標系下的彈丸質(zhì)心運動方程[3]可表示為

(1)

(2)

根據(jù)前文坐標系和角度的定義可知，

(3)

(4)

式中：Fyb、Fzb分別表示合外力在彈道坐標系內(nèi)Oyb軸和Ozb軸上的投影；ωηf、ωζf分別表示彈丸擺動角速度在第一彈軸系內(nèi)Oηf軸和Oζf軸上的投影。

結(jié)合彈軸擺動的動力學方程[3]，可得旋轉(zhuǎn)穩(wěn)定彈丸非線性角運動方程為

(5)

式中：ωξ表示彈丸轉(zhuǎn)動角速度在第一彈軸系內(nèi)Oξ軸上的投影；A和C分別表示彈丸的赤道轉(zhuǎn)動慣量和極轉(zhuǎn)動慣量；Mηf、Mζf分別表示彈丸所受合外力矩在第一彈軸系內(nèi)Oηf軸和Oζf軸上的投影。

1.3 旋轉(zhuǎn)穩(wěn)定彈丸氣動力模型

這里只介紹非線性角運動方程中需要用到的氣動力模型。

1.3.1 升力

升力在彈道坐標系內(nèi)Oyb軸和Ozb軸上的投影表達式為

(6)

式中：ρ為空氣密度；S為參考面積；Cl為彈丸升力系數(shù)。

1.3.2 靜力矩

靜力矩在第一彈軸系內(nèi)Oηf軸和Oζf軸上的投影表達式為

(7)

式中：l為參考長度；mz為靜力矩系數(shù)。

1.3.3 赤道阻尼力矩

赤道阻尼力矩在第一彈軸系內(nèi)Oηf軸和Oζf軸上的投影表達式為

(8)

式中：d為彈徑；mzz為赤道阻尼力矩系數(shù)。

1.3.4 馬格努斯力矩

馬格努斯力矩的非線性現(xiàn)象在很多實驗中被發(fā)現(xiàn)。文獻[2]研究表明，非線性馬格努斯力矩是影響旋轉(zhuǎn)彈穩(wěn)定性的主要因素之一。因此本文考慮馬格努斯力矩的非線性項，并建立模型如下：

(9)

式中：my為馬格努斯力矩系數(shù)，

(10)

my0、my2分別表示馬格努斯力矩系數(shù)一次項和三次項系數(shù)，δe為彈箭總攻角，cosδe=cosδacosδd.

由于重力只影響系統(tǒng)的平衡點位置，且通過計算可知這種影響非常小。另外非線性系統(tǒng)的平衡點可以通過坐標變換移到原點，因此重力影響可忽略。綜上所述，可得

(11)

2 平衡點性質(zhì)分析

對于(5)式所表示的非線性系統(tǒng)，令方程組等號右邊為0，可得原點[0 rad 0 rad 0 rad/s 0 rad/s]T是其一個孤立平衡點。根據(jù)三角函數(shù)的周期性可知，系統(tǒng)有9個平衡點，可表示為

[k1π rad,k2π rad,0 rad/s,0 rad/s]T,

(12)

式中：k1、k2可取0和±1.

系統(tǒng)在平衡點處的穩(wěn)定性可在平衡點處對方程線性化，再通過分析其雅克比矩陣(見(13)式)的性質(zhì)來判定:

(13)

式中：Q為系統(tǒng)的雅克比矩陣；fi(i=1,2,3,4)為系統(tǒng)第i個狀態(tài)方程；xi為系統(tǒng)的第i個狀態(tài)量。

因此系統(tǒng)在平衡點處的雅克比矩陣可表示為

(14)

通過計算和仿真可知，系統(tǒng)的平衡點性質(zhì)和馬格努斯力矩三次項系數(shù)my2有關。當my2=0時，平衡點的穩(wěn)定性如圖1所示，穩(wěn)定的平衡點和不穩(wěn)定的平衡點相互間隔。此時，所有穩(wěn)定的平衡點都是穩(wěn)定的焦點，相互間是無差別的。根據(jù)對稱性可知，原點的吸引域為δa,δd∈[-π/2 rad,π/2 rad]。如果初始δa和δd在該范圍內(nèi)，且角速度不是很大，則系統(tǒng)最終收斂于原點，如圖2(a)所示。如果初始角速度足夠大，則會使兩個攻角超出吸引域的范圍，被其他穩(wěn)定的平衡點所吸引，如圖2(b)所示。對于my2=0的情況，平衡點的吸引域不受系統(tǒng)狀態(tài)參數(shù)的影響，無論對于什么發(fā)射條件，平衡點的吸引域大小和形狀都不會發(fā)生改變。

圖1 平衡點穩(wěn)定性示意圖

圖2 不同初值下系統(tǒng)的收斂情況

當my2≠0時，平衡點的穩(wěn)定性和穩(wěn)定平衡點的吸引域大小受系統(tǒng)參數(shù)影響。通過仿真計算可知，此時平衡點的穩(wěn)定性分為以下兩種情況：

情況1系統(tǒng)有5個穩(wěn)定的平衡點，分別為原點[0 rad 0 rad 0 rad/s 0 rad/s]T、[±π rad 0 rad 0 rad/s 0 rad/s]T和[0 rad ±π rad 0 rad/s 0 rad/s]T，其余平衡點均不穩(wěn)定。

情況2系統(tǒng)只有原點一個穩(wěn)定的平衡點，其余平衡點均不穩(wěn)定。

情況1與my2=0的情況相似，原點與其周圍4個相鄰的平衡點性質(zhì)一致，這樣原點的吸引域亦為δa,δd∈[-π/2 rad,π/2 rad]。對于情況2，如果此時系統(tǒng)的初始點在不穩(wěn)定平衡點附近，則系統(tǒng)就會不穩(wěn)定。如果初始點在原點附近，則系統(tǒng)將保持穩(wěn)定。這表明原點周圍存在一個吸引域，且該吸引域的邊界是一個不穩(wěn)定的極限環(huán)。

根據(jù)非線性穩(wěn)定性理論，不穩(wěn)定的極限環(huán)內(nèi)包含一個穩(wěn)定的平衡點，屬于次臨界Hopf分岔的一個狀態(tài)。因此，對于這一類平衡點吸引域邊界的計算就可以轉(zhuǎn)換為次臨界Hopf分岔狀態(tài)下極限環(huán)半徑的求解?？梢匀芜x一個系統(tǒng)參數(shù)作為分岔參數(shù)，計算分岔點。在分岔點處，對系統(tǒng)方程使用中心流形定理進行降維處理，然后對降維后的系統(tǒng)使用次臨界Hopf分析，就可獲得極限環(huán)的半徑。這樣也就求得了該平衡點的吸引域邊界。

3 平衡點吸引域邊界計算

3.1 構造分岔條件

中心流形定理規(guī)定，若系統(tǒng)雅克比矩陣Q的一部分特征值實部為0，其余特征值具有負實部，此時無法使用線性化來確定原點的穩(wěn)定性，而系統(tǒng)可以轉(zhuǎn)化為1個階數(shù)與特征值為負實部的個數(shù)一致的降階系統(tǒng)[14]。本文通過構造分岔條件，從而達到滿足使用中心定理的條件。

任選系統(tǒng)一個參數(shù)μ為分岔參數(shù)，然后求出帶參數(shù)μ的矩陣Q的特征多項式

f(λ,μ)=a4λ4+a3λ3+a2λ2+a1λ+a0,

(15)

式中：λ為特征值；a0、a1、a2、a3、a4均為含有參數(shù)μ的函數(shù)。根據(jù)4階系統(tǒng)的霍爾維茨判據(jù)可知，要想4階系統(tǒng)穩(wěn)定，必須使得(15)式中的所有系數(shù)大于0，且其霍爾維茨行列式滿足(16)式條件：

(16)

式中：Δ2、Δ3分別為矩陣Q的2階和3階主子式。

因此只需尋找參數(shù)值μ0，使得(16)式等于0，且(15)式中的所有系數(shù)仍大于0，即滿足中心定理的要求。此時參數(shù)值μ0也為次臨界Hopf分岔中的分岔點。

3.2 非線性角運動方程極限環(huán)的判定和計算

在分岔點μ0處，系統(tǒng)滿足使用中心流形定理的要求，(15)式有一對實部為0的特征根和一對實部為負的特征根。因此，系統(tǒng)(5)式存在二維中心流形，可轉(zhuǎn)換為二維系統(tǒng)。最終得到的中心流形上流的約化方程可表示為

(17)

式中：y1、y2為原系統(tǒng)狀態(tài)變量經(jīng)非奇異線性變換后的狀態(tài)量。

(17)式寫成如下形式：

(18)

式中：α(μ)、β(μ)分別為零點鄰域的線性近似系統(tǒng)共軛特征值的實部和虛部；g1(y1,y2,μ)、g2(y1,y2,μ)分別表示f1(y1,y2)和f2(y1,y2)經(jīng)非奇異線性變換后的3階小量。

根據(jù)規(guī)范性理論[15]，寫出(18)式的3階規(guī)范形為

(19)

式中：

(20)

μ0為分岔點；ω為β(μ)在分岔點μ0處的取值；e為?β(μ)/?μ在分岔點μ0處的取值；a為中間變量，

(21)

b為規(guī)范變換的中間變量，由于b對平衡點的分岔特性無影響且推導過程復雜，本文參照文獻[15]的處理方法，不推導其具體表達形式。

當a≠0、c≠0時，系統(tǒng)在μ=μ0處出現(xiàn)Hopf分岔，根據(jù)a和c的符號就可以判定分岔特性。因為本文構造的是次臨界Hopf分岔，故所得結(jié)果滿足a>0、c>0，原點對μ<μ0漸進穩(wěn)定，且此時存在不穩(wěn)定的極限環(huán)。該極限環(huán)就是漸進穩(wěn)定平衡點吸引域的邊界，其半徑可由(22)式求得

(22)

此時系統(tǒng)的分岔圖如圖3所示。

圖3 系統(tǒng)分岔圖

4 各參數(shù)對角運動影響的計算分析

3.2節(jié)給出了原點平衡點吸引域的計算方法，現(xiàn)以某型155 mm榴彈為例，通過計算分析各參數(shù)對彈丸穩(wěn)定性的影響，從而解釋一些特殊的彈丸運動現(xiàn)象。計算所采用的彈丸主要特性參數(shù)和氣動力數(shù)據(jù)分別如表1和表2所示。

表1 彈丸仿真參數(shù)

表2 彈丸氣動參數(shù)

4.1 馬格努斯力矩三次項系數(shù)的影響

考慮彈丸初速為510 m/s時，不同馬格努斯力矩三次項系數(shù)對原點平衡點吸引域邊界的影響，計算結(jié)果如圖4所示。

圖4 不同my2值的不穩(wěn)定極限環(huán)

由圖4可知，my2的絕對值越大，原點的吸引域就越小，因此在進行彈丸氣動力設計時，應盡量保證my2的絕對值較小。圖4中my2=0.02時，平衡點穩(wěn)定性屬于情況1，通過第2節(jié)分析可知，其吸引域應為δa,δd∈[-π/2 rad,π/2 rad]。此時a<0、c>0，平衡點周圍無不穩(wěn)定的極限環(huán)。由圖4(b)可以看出，當my2為正且平衡點穩(wěn)定性屬于情況2時，吸引域面積隨著my2的增大快速減小，因此在設計時也應避免這種情況的發(fā)生。

圖5給出了my2=0.2時狀態(tài)初值分別在吸引域內(nèi)和吸引域外的仿真結(jié)果。由圖5可知，當狀態(tài)初值在吸引域內(nèi)部時，系統(tǒng)會正常收斂于原點。當狀態(tài)初值在吸引域外部時系統(tǒng)會發(fā)散，出現(xiàn)運動失穩(wěn)現(xiàn)象。

圖5 初值在吸引域內(nèi)外的仿真結(jié)果

4.2 初速的影響

考慮my2=0.2時，不同初速對原點吸引域邊界的影響，計算結(jié)果如圖6所示。

圖6 不同初速值的不穩(wěn)定極限環(huán)

從圖6中可以看出，當初速從小到大變化時，原點的吸引域面積先減小后增大。這就能解釋為什么“彈道炸”和美國M549彈丸發(fā)射失穩(wěn)事故總是發(fā)生在某一特定的裝藥號下。在這些裝藥號下，原點的吸引域相對于其他初速是最小的，因此是最有可能出現(xiàn)運動的失穩(wěn)。

4.3 空氣密度的影響

考慮my2=2.0、速度為510 m/s時，不同空氣密度對原點平衡點吸引域邊界的影響，計算結(jié)果如圖7所示。

圖7 不同密度值的不穩(wěn)定極限環(huán)

從圖7中可以看出，當密度逐漸增大時，原點的吸引域面積在逐漸減小，表明大的空氣密度不利于旋轉(zhuǎn)彈的穩(wěn)定。這一結(jié)論與文獻[7]一致，但二者表述方式不同。文獻[7]認為，旋轉(zhuǎn)穩(wěn)定彈丸如果在平原(空氣密度大)使用時滿足動態(tài)穩(wěn)定條件，則在高原(空氣密度小)使用時更容易滿足動態(tài)穩(wěn)定性條件。這是一種線性穩(wěn)定性條件，只與彈丸特征參數(shù)和氣動力系數(shù)有關，與彈丸角運動狀態(tài)量的初始值無關。本文通過計算獲得密度與原點吸引域大小的關系，從而將彈丸的穩(wěn)定性與彈丸發(fā)射條件聯(lián)系起來。只要彈丸的角運動不超過吸引域范圍，則彈丸可保持穩(wěn)定飛行。

空氣密度對原點吸引域大小的影響規(guī)律也可用于解釋一些旋轉(zhuǎn)彈丸運動不穩(wěn)定的現(xiàn)象。美國M549彈丸發(fā)射失穩(wěn)全部發(fā)生于寒冷的冬季，是因為冷空氣的密度大，原點吸引域變小，彈丸的起始擾動很容易使得彈丸攻角超出吸引域之外，從而發(fā)生運動失穩(wěn)現(xiàn)象。“彈道炸”事故多發(fā)生于小射角下，是因為小射角發(fā)射的彈丸彈道平直，高度變化小，從而密度始終保持在較大值，因此更容易發(fā)生“彈道炸”事故。

4.4 特殊發(fā)射條件下的穩(wěn)定性分析

移動發(fā)射平臺(如艦船、坦克等)在行進間發(fā)射彈丸時，由于發(fā)射平臺具有一定大小的運動速度，出炮口的彈丸速度方向與身管方向不一致。另一方面，彈丸受彈帶的約束作用，彈軸方向基本和身管方向保持一致。因此，高速行駛的艦船兩側(cè)發(fā)射的彈丸方向攻角δd大小相等、方向相反，如圖8所示。根據(jù)1.1節(jié)的定義可知，彈軸在速度軸右側(cè)時方向攻角為正，故左側(cè)彈丸攻角δl為負，右側(cè)彈丸攻角δr為正。

圖8 艦炮發(fā)射彈丸方向攻角示意圖

彈丸在出炮口時，由于重力的傾離作用，使彈軸具有繞彈帶向下擺動的角速度，即ωζf為負。對于艦炮兩側(cè)發(fā)射的彈丸，該角速度大小相等、方向相同，均指向Oζf軸負向。

如果彈丸原點的吸引域比較小，而起始攻角又相對較大，則對于不同初始條件，彈丸的穩(wěn)定性也會出現(xiàn)不同情況。圖9顯示了不同起始條件彈丸角運動的相圖。仿真條件為：my2=2.0、初速510 m/s、密度1.2 kg/m3、左舷發(fā)射初始狀態(tài)量為[0° -6° 0 rad/s -2 rad/s]T、右舷發(fā)射初始狀態(tài)量為[0° 6° 0 rad/s -2 rad/s]T。

圖9 艦炮左右舷發(fā)射彈丸穩(wěn)定性示意圖

由圖9可知，對于行進間艦船，左舷發(fā)射的旋轉(zhuǎn)彈丸不穩(wěn)定，右舷發(fā)射的旋轉(zhuǎn)彈丸穩(wěn)定。圖9中彈丸的起始攻角都在原點吸引域以內(nèi)，且起始角速度大小相同，彈丸角運動的穩(wěn)定性卻不同。當起始方向攻角δd和起始彈軸擺動角速度ωζf同向時，δd開始會增大。當起始方向攻角δd和起始彈軸擺動角速度ωζf反向時，δd開始會減小。因此，當起始方向攻角比較大時，同向的ωζf就有可能使δd超出原點的吸引域，最終導致彈丸不穩(wěn)定。

5 結(jié)論

本文推導了旋轉(zhuǎn)穩(wěn)定彈丸的非線性角運動模型，針對周期性平衡點分類分析了原點平衡點的吸引域。針對原點為唯一穩(wěn)定平衡點的情況，通過計算原點處次臨界Hopf分岔極限環(huán)半徑得到其吸引域，并計算了多種條件對吸引域的影響。得出以下主要結(jié)論：

1)馬格努斯力矩的三次項系數(shù)對原點吸引域大小影響最大，非線性馬格努斯力矩是影響旋轉(zhuǎn)穩(wěn)定彈丸非線性穩(wěn)定性的主要因素之一。

2)在超音速范圍內(nèi)，隨著彈丸初速的增大，原點吸引域的面積先減小后增大，因此在發(fā)射時存在一個最差初速，在進行裝藥號設計時要予以考慮，盡量規(guī)避該初速。

3)空氣密度越大，彈丸的原點吸引域面積越小，解釋了在冬季小射角射擊時的彈丸失穩(wěn)現(xiàn)象。

4)分析了移動發(fā)射平臺右側(cè)發(fā)射的彈丸穩(wěn)定，而左側(cè)發(fā)射彈丸不穩(wěn)定的原因：兩側(cè)發(fā)射彈丸的方向攻角方向相反，在相同彈軸擺動角速度的作用下，右側(cè)方向攻角減小，收斂于原點；左側(cè)方向攻角會增大超出原點吸引域，造成飛行不穩(wěn)定。