旋轉穩(wěn)定彈丸非線性角運動吸引域計算方法

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

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

0 引言

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

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

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

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

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

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

1.1 坐標系和角度的定義

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

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

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

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

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

(1)

(2)

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

(3)

(4)

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

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

(5)

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

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

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

1.3.1 升力

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

(6)

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

1.3.2 靜力矩

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

(7)

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

1.3.3 赤道阻尼力矩

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

(8)

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

1.3.4 馬格努斯力矩

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

(9)

式中：my為馬格努斯力矩系數，

(10)

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

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

(11)

2 平衡點性質分析

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

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

(12)

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

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

(13)

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

因此系統在平衡點處的雅克比矩陣可表示為

(14)

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

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

圖2 不同初值下系統的收斂情況

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

情況1系統有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系統只有原點一個穩(wěn)定的平衡點，其余平衡點均不穩(wěn)定。

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

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

3 平衡點吸引域邊界計算

3.1 構造分岔條件

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

任選系統一個參數μ為分岔參數，然后求出帶參數μ的矩陣Q的特征多項式

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

(15)

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

(16)

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

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

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

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

(17)

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

(17)式寫成如下形式：

(18)

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

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

(19)

式中：

(20)

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

(21)

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

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

(22)

此時系統的分岔圖如圖3所示。

圖3 系統分岔圖

4 各參數對角運動影響的計算分析

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

表1 彈丸仿真參數

表2 彈丸氣動參數

4.1 馬格努斯力矩三次項系數的影響

考慮彈丸初速為510 m/s時，不同馬格努斯力矩三次項系數對原點平衡點吸引域邊界的影響，計算結果如圖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)初值分別在吸引域內和吸引域外的仿真結果。由圖5可知，當狀態(tài)初值在吸引域內部時，系統會正常收斂于原點。當狀態(tài)初值在吸引域外部時系統會發(fā)散，出現運動失穩(wěn)現象。

圖5 初值在吸引域內外的仿真結果

4.2 初速的影響

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

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

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

4.3 空氣密度的影響

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

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

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

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

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

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

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

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

如果彈丸原點的吸引域比較小，而起始攻角又相對較大，則對于不同初始條件，彈丸的穩(wě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ā)射的旋轉彈丸不穩(wěn)定，右舷發(fā)射的旋轉彈丸穩(wěn)定。圖9中彈丸的起始攻角都在原點吸引域以內，且起始角速度大小相同，彈丸角運動的穩(wěn)定性卻不同。當起始方向攻角δd和起始彈軸擺動角速度ωζf同向時，δd開始會增大。當起始方向攻角δd和起始彈軸擺動角速度ωζf反向時，δd開始會減小。因此，當起始方向攻角比較大時，同向的ωζf就有可能使δd超出原點的吸引域，最終導致彈丸不穩(wěn)定。

5 結論

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

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

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

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

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