旋轉穩定彈丸非線性角運動吸引域計算方法

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

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

0 引言

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

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

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

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

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

1 旋轉穩定彈丸動力學建模

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 旋轉穩定彈丸非線性角運動模型

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

(1)

(2)

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

(3)

(4)

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

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

(5)

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

1.3 旋轉穩定彈丸氣動力模型

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

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 馬格努斯力矩

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

(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.

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

(13)

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

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

(14)

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

圖1 平衡點穩定性示意圖

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

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

情況1系統有5個穩定的平衡點，分別為原點[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，其余平衡點均不穩定。

情況2系統只有原點一個穩定的平衡點，其余平衡點均不穩定。

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

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

3 平衡點吸引域邊界計算

3.1 構造分岔條件

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

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

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

(15)

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

(16)

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

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

3.2 非線性角運動方程極限環的判定和計算

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

(17)

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

(17)式寫成如下形式：

(18)

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

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

(19)

式中：

(20)

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

(21)

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

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

(22)

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

圖3 系統分岔圖

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

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

表1 彈丸仿真參數

表2 彈丸氣動參數

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

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

圖4 不同my2值的不穩定極限環

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

圖5給出了my2=0.2時狀態初值分別在吸引域內和吸引域外的仿真結果。由圖5可知，當狀態初值在吸引域內部時，系統會正常收斂于原點。當狀態初值在吸引域外部時系統會發散，出現運動失穩現象。

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

4.2 初速的影響

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

圖6 不同初速值的不穩定極限環

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

4.3 空氣密度的影響

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

圖7 不同密度值的不穩定極限環

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

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

4.4 特殊發射條件下的穩定性分析

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

圖8 艦炮發射彈丸方向攻角示意圖

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

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

圖9 艦炮左右舷發射彈丸穩定性示意圖

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

5 結論

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

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

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

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

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