基于歐拉法的自尋的炮射導彈彈道建模與仿真

白 卓，趙河明，楊晉偉

（1.中北大學機電工程學院，太原030051；2.淮海工業(yè)集團，山西長治046012）

常規(guī)武器發(fā)射平臺在使用和維護不變的情況下，遇到要求精確打擊時，可直接使用制導彈藥。這種擴展使常規(guī)武器兼?zhèn)淞顺Ｒ?guī)武器的傳統(tǒng)功能及過去只有制導武器才具有的精確打擊功能。從20 世紀的幾次局部戰(zhàn)爭來看，精確制導武器在現(xiàn)代戰(zhàn)爭中起著決定性的作用。據有關資料統(tǒng)計，過去平均使用250 發(fā)155 mm 的常規(guī)彈藥，只可能擊毀一輛坦克，現(xiàn)在使用精確制導技術的炮彈僅需1～2 發(fā)即可，其效能提高了125～259倍?，F(xiàn)在世界各國都在爭相研制的炮射導彈就是一種精確制導武器[1-7]。裝有精確制導系統(tǒng)的炮射導彈的出現(xiàn)使現(xiàn)代火炮進入了一個新的時代，實現(xiàn)了彈炮集合，由無控向有控的轉變[8]。相對于機載導彈，炮射導彈目前還有很大的發(fā)展?jié)摿?，特別是炮射自尋的防空導彈技術，使坦克擁有抗衡武裝直升機，甚至攻擊機的“撒手锏”[9]。

大多數(shù)炮射導彈采用直瞄攻擊模式，不僅耗時長，而且精確度與毀傷無法達到預想效果。目前，西方國家制導彈藥多采用制導體制，自尋的模式代表著制導彈藥的發(fā)展方向。

1 炮射導彈系統(tǒng)數(shù)學模型

由理論力學可知，導彈運動在每一瞬時，把導彈當作一個質量不變的，在氣動力、推力、操縱力等作用下運動的剛體來處理[9]。

通過這些假定，來研究導彈質心的移動。并且，在利用質心運動的動量定理時，應將彈道坐標系作為動坐標系，這樣會使方程式簡單明了，更便于編寫程序，求解方程組。

1.1 導彈的受力分析

導彈在地球周圍大氣里飛行，會受到重力、空氣動力和發(fā)動機推力的共同作用。在導彈上的重力為地心引力和離心慣性力的矢量和。

式（1）、（2）中：cx、cy、cz分別為阻力系數(shù)、升力系數(shù)和側向力系數(shù)；mx、my、mz分別稱為滾動力矩系數(shù)、偏航力矩系數(shù)和俯仰力矩系數(shù)；L為特征長度；S為特征面積。

發(fā)動機推力是由發(fā)動機內的燃氣流以高速噴射出而產生的反作用力。

1.2 導彈數(shù)學模型建立

本文的彈道設計為方案飛行彈道，末端采用比例導引法制導，實現(xiàn)炮射導彈的“自尋的”。導彈的飛行路線設計：炮口發(fā)射→無控階段→程序爬升階段→平飛階段→比例導引段毀傷目標。圖1為導彈和目標的運動學幾何關系。

圖1 導彈和目標的運動學幾何關系Fig.1 The kinematic geometry relationshipbetween missile and target

1）無控升段。在無控彈道段為自由飛行，忽略其他干擾，只考慮重力和氣動力的作用。

2）程序爬升階段。對于程序爬升段，目的是提升彈道高度，增加射程和為實現(xiàn)大落腳準備。對于爬升段的方案飛行，考慮到技術上的可行性，實現(xiàn)的難易程度以及精確性，本文采用給定俯仰角變化規(guī)律??(t)的飛行方案。這段時間內，發(fā)動機開始工作，相比于無控段的彈道方程組，爬升段要引入一個新的變量mc（質量秒流量）。

3）平飛階段。導彈在平飛段不會進行大的機動飛行（速度、高度恒定），并且在飛行期間彈道傾角不會改變。對于平飛階段，常采用超低空飛行或大高度飛行，以避開敵人雷達偵查和防空火力的射擊。

式（3）中：v為導彈的速度；P為發(fā)動機推力；X為軸向力；g為重力加速度；θ為彈道傾角；ωz為俯仰角角速度；Mz為靜俯仰力矩；Jz為赤道轉動慣量；-mc為質量流量；?為俯仰角；α為平衡攻角。

導彈—目標相對位置：

導彈—目標相對速度：

導彈—目標相對距離：

導彈—目標視線角速度：

導彈—目標相對距離變化率：

2 彈道方程數(shù)值解法

彈道方程組是一階變系數(shù)聯(lián)立方程組，只能用數(shù)值方法求得數(shù)值解，僅在一些特定條件下經過適當?shù)暮喕拍芮蟮媒平馕鼋狻?/p>

本節(jié)通過對比龍格—庫塔法（Runge—Kutta）、阿當姆斯法以及歐拉法，選出合適的方法求解彈道微分方程[13-16]。

2.1 龍格—庫塔法

龍格—庫塔法實質上是以函數(shù)y(x)的臺勞級數(shù)為基礎的一種改進方法。最常用的是4階龍格—庫塔法，若已知在點n處的值(tn,y1n,y2n,…,ymn) ，則求點n+1處的函數(shù)值的公式可由龍格—庫塔公式推導出：

式（9）中：

2.2 阿當姆斯預報—校正法

2.3 歐拉法

3 自尋的炮射導彈彈道仿真

3.1 仿真方法

本文采用歐拉法來求解彈道方程組，在初始條件給定的情況下，程序可以得到以時間t為自變量的彈道上的任意諸元。對于剛體彈道方程，時間的步長h必須小于0.005 s。否則，計算發(fā)散。故在使用歐拉法求解彈道方程組時，選取0.001作為積分步長，既可便捷的求解彈道方程，又可滿足誤差要求。比例導引方法的比例系數(shù)K，應選擇在1～∞的范圍內，通常在2～6 的范圍內選擇。本文選擇K=4。

采用Matlab 編寫程序來實現(xiàn)自尋的炮射導彈的各段彈道，程序以時間t作為全局變量，將計算出的數(shù)值存儲在數(shù)組中，便于以后作圖[17-20]。

3.2 仿真結果及分析

仿真初始條件為：炮射導彈以初速度300 m/s 發(fā)射，初始攻角α=0°，初始俯仰角?=18°，初始彈道的傾角θ=18°，得到仿真結果如圖2所示。

圖2 彈道曲線Fig.2 Ballistic curve

圖2為自尋的炮射導彈的彈道曲線。從圖2可以得出，導彈以初速300 m/s 發(fā)射，在發(fā)動機不工作的情況下爬升一段距離。此時，彈道傾角逐漸減?。浑S后，發(fā)動機開始工作，導彈在發(fā)動機的推動下爬升，爬升過程中隨著彈道傾角的增加，爬升率下降；當t=22 s時，導彈爬升到最高點；然后，經歷一段平飛，導彈迅速調整姿態(tài)，彈道傾角在平飛段為0，平飛段的目的為使導彈飛到目標上空，增加導彈射程；當彈目距離達到預定值之后，導彈加速下落，比例導引法擊毀目標。

平衡攻角和俯仰角與彈道傾角之間存在著簡單的幾何關系：α=?-θ。

如圖3 所示，在導彈發(fā)射時，平衡攻角α=0°，平衡攻角α一直保持為0°，導彈保持平衡飛行，直至彈道轉為比例導引段，導彈姿態(tài)調整，進入末端制導，實際攻角開始減小，其收斂于平衡攻角，導彈飛行穩(wěn)定。顯然，圖示攻角一直保持為0，僅為理想狀態(tài)下攻角曲線，實際的攻角曲線應如圖4所示。

圖4 所示為某無控彈道的攻角的仿真曲線，舵偏角δz=0°，與此對應的平衡攻角α也為0。導彈在飛行過程中攻角的變化會產生嚴重的振蕩，因為彈體阻尼力矩的存在，攻角的變化會收斂于平衡攻角。

圖3 平衡攻角-時間曲線Fig.3 Angle of attack-time curve

圖4 平衡攻角-時間曲線Fig.4 Angle of attack-time curve

由圖5 速度時間變化的曲線可以看出，速度曲線近似為3段折線。導彈經發(fā)射之后，在0～30 s，導彈在較高的切向加速度加速飛行；在30 s 之后，導彈到達預定飛行高度，速度保持在625 m/s 左右，開始以勻速飛行；150 s 之后，導彈進入末端制導段，速度開始減小，飛向并擊毀目標。

依據圖6 的彈道傾角隨時間的變化曲線可以看出，在無控階段，彈道傾角逐漸減小，直到發(fā)動機開始工作，彈道的傾角開始增大；導彈轉至平飛段后，彈道傾角急速減小到0，并在平飛段一直保持在0；在比例導引階段，導彈加速向目標飛去，并且在此過程中，彈道傾角不斷減小。

圖5 速度-時間曲線Fig.5 Velocity-time curve

圖6 彈道傾角-時間曲線Fig.6 Ballistic inclination-time curve

圖7 俯仰角-時間曲線Fig.7 Pitch angle-time curve

4 結論

根據以上一系列圖示的仿真結果，炮射導彈可以實現(xiàn)自尋的模式。本文在前人提出的“自尋的”炮射導彈的基礎之上，設計了上文所示的可以實現(xiàn)頂部攻擊的、對靜止目標或慢速移動目標實現(xiàn)高毀傷打擊的炮射導彈彈道模型。

從仿真的結果來看，當初速等參數(shù)一定時，就可以確定一條彈道曲線。所以，只要改變模型中的相應參數(shù)，導彈飛行時間、射程、彈道高度、落點速度都會做出相應的改變。本文所提到的研究方法、思路和分析方法可以為未來實現(xiàn)炮射導彈智能化提供借鑒，具有一定的實際意義。