車載導彈推力偏心初始擾動影響研究*

鐘洲，孟令濤，曾偉，曾魯，王晨曦

(中國運載火箭技術研究院，北京 100076)

0 引言

初始擾動是影響導彈發射精度的主要因素之一，是由發射過程中的有關干擾因素引起的，這些干擾因素涉及到武器系統的各個方面，一般都屬于隨機性質，包括質量偏心、動不平衡、運載發射環境等多種激勵。分析隨機激勵下導彈姿態的分布有助于理解隨機激勵與導彈初始擾動的內在關系，提高研究人員對隨機激勵下的發射動態響應規律的認識，有助于發射系統的優化設計[1-2]。

國外學者[3]通過建立含結構特性干擾因素的動力學和運動學方程開展了發射動力學初始擾動分析。國內陳陣[4]等通過建立半約束期火箭彈的運動微分方程，研究了質量偏心和動不平衡等隨機激勵對火箭初始擾動的影響；李翔[5]等通過設計虛擬正交試驗和多因素極差分析研究了艦船搖擺對導彈初始擾動的影響；康甜[6]等通過建立引入干擾因素的偏差模型分析了推力偏心、導軌不平度對初始擾動的影響。發動機由于裝藥燃燒特性和噴管垂直度等原因，其產生的推力存在偏心現象，不僅會對飛行階段的平穩控制產生重要隱患[7]，亦會對導彈起飛階段的姿態產生一定影響。因此，本文主要研究發動機推力偏心對車載導彈垂直發射時的初始擾動影響，利用Lagrange方程和修正的Craig-Bampton法建立了車載導彈剛柔耦合垂直發射動力學模型；通過參數化腳本語言(.cmd)[8]和C++編譯程序，進行了參數化動力學仿真，獲得大量不同推力偏心情況下的導彈初始擾動，并進一步運用直方圖法[9]和極大似然估計法對樣本數據進行了數據統計分析，研究了導彈姿態的分布規律，得到了有益的初步結論，為隨機激勵的進一步研究提供了一定的基礎。

1 基于修正的Craig-Bampton法多柔體動力學理論

1.1 修正的Craig-Bampton

結構動力學中傳統的Craig-Bampton(簡稱C-B)方法[10]是基于彈性體沒有大范圍剛體運動的假設，分析的是線性動力學問題。子結構裝配時，嵌入到C-B約束模態中包含6個剛體自由度，并且約束模態的本質是靜力縮聚情況下計算出的模態，其模態和頻率不能相互對應。而在多體系統動力學中，已經定義了構件的剛體位移，同時構件在慣性坐標系中有大范圍運動，屬于非線性動力學問題，因此必須對傳統Craig-Bampton子結構法進行適當修正[11]，使其能滿足多體系統動力學計算的需要。

易知，在C-B模態坐標下，子結構的結構動力學方程為

(1)

式中：Mi,Ci,Ki,Ri分別為C-B模態基下的子結構質量、阻尼、剛度和外力矩陣。

求解式(6)所對應無阻尼振動方程的特征值和特征向量：

(2)

(3)

因此，子結構的物理自由度可以表示為

(4)

利用上述正交轉換修正后，使得對界面坐標依次固定時產生的模態，被無約束的模態近似代替；約束模態被邊界特征向量所取代。因此，可以除去零頻率對應的6個剛體運動模態，并且所有模態和頻率能相互對應。

1.2 多體系統剛柔耦合動力學方程

如果物體坐標系的位置，使用它在慣性參考系中的笛卡爾坐標X=(x,y,z)T和反映剛體方位的歐拉角ψ=(ψ,θ,φ)T來表示，并且模態坐標用式q=(q1,…,qM)T表示(M為模態坐標數)，那么柔性體的廣義坐標可以表示為ξ=(X,ψ,q)T。

運用拉格朗日方程可以建立柔性體的運動微分方程：

(5)

式中：K和D分別為模態剛度矩陣和模態阻尼矩陣；Kξ和Dξ分別為物體內部由于彈性變形和阻尼引起的廣義力；fg為廣義重力；λ對應于約束的拉格朗日乘子；Q為對應于外力的廣義力。

將柔性體的計算結果和多剛體理論方法相結合可以進一步得到多體系統的剛柔耦合動力學方程。通過文獻[12-14]可以推導出使用拉格朗日乘子得到的第i個柔體或剛體的方程形式為

(6)

式中：K為動能的表達式;Qi為廣義力，包括單元彈性變形和外加載荷(包括單元間載荷)引起的廣義力。將式(11)與約束方程C(q,t)=0聯立，即構成剛柔耦合多體動力學方程。

2 車載導彈剛柔耦合動力學模型

2.1 武器系統組成分析與假設

車載導彈垂直發射系統是一個復雜的多體機械系統，發射時導彈垂直裸立在發射臺支撐盤上，發射載荷由發射臺傳遞至車體，其動力學模型為考慮接觸碰撞效應的多個結構件(剛體或柔體)通過各種運動副、等效彈簧、阻尼器連接組成的系統。同時本文研究的對象采用半剛性支撐方式，即發射時由輪胎和調平支腿共同承擔發射載荷，可見其垂直發射過程力學環境非常復雜，相當于一個多因素耦合系統，具有高非線性特點。因此，建立能反映發射動態特性的柔性多體數學模型是進行初始擾動分析的關鍵因素。

根據車載導彈武器系統的子部件功能關系，對整車系統進行合理處理，簡化后的模型主要分為車輪、調平支腿、懸掛系統、車體、上裝設備、導流臺、發射臺、起豎油缸和導彈。其中，導彈無約束豎立在發射臺支撐盤上；發射臺通過旋轉鉸與車體和起豎油缸連接；導流臺、車輪和調平支腿均與地面接觸；上裝設備緊固在車體上；車體與車輪通過懸掛系統連接。系統主要部分具體拓撲連接如圖1所示。

為了描述武器系統垂直發射時的動力學特性，特作如下假設：

(1) 不考慮連接鉸的間隙，均設為理想約束；

(2) 懸掛系統簡化成旋轉或平移的彈性元件，利用線性和非線性曲線描述其剛度和阻尼特性(仿真中利用Vforce 力元模擬)；

(3) 考慮地面為理想的剛體；

(4) 不考慮發射過程中導彈的變質量特性。

2.2 柔性體建模

整車模型中關鍵結構件的柔性效應會對導彈姿態產生重要影響，需要利用相關技術手段完成柔性體建模。考慮到系統各結構特性對發射過程的影響程度，將發射臺作為柔性體。Adams/Flex是采用修正的Craig-Bampton法來描述彈性體的變形，同時它給Adams和Abaqus提供了雙向數據交換接口，本文利用有限元軟件Abaqus對發射臺實體模型進行網格劃分，定義好附著點(施加力與約束副的外部節點)和構件間實際連接處節點的關系，通過宏命令得到模態中性文件(.mnf)，導入Adams生成柔性體提高仿真的精確度。發射臺柔性體有限元模型材料特性密度為7.8×103kg/m3，彈性模量為210 GPa，泊松比為0.3，共有6個附著點和324 578個節點數。表1列出了其仿真主要模態頻率。

圖1 拓撲連接圖Fig.1 Topological connection graph

模態階數789101112模態頻率/Hz35.943.252.562.880.7102.8

2.3 激勵

發射過程中導彈和發射臺的接觸力用contact(solid body to flex)工具定義，其本質為：接觸碰撞采用等效彈簧阻尼模型，即彈簧接觸力根據Hertz接觸理論計算，同時用阻尼器模擬接觸過程中的能量損失；間隙處的切向摩擦力采用coulomb摩擦模型計算。Contact工具中的具體參數值與接觸材料有關，可根據試驗數據和軟件提供的經驗數據設定：接觸剛度為3.8×104N/mm；阻尼系數為28 N·s/mm；力貢獻指數為1.5；最大穿透深度0.1 mm；靜摩擦系數為0.3；動摩擦系數為0.1。

導彈推力由發動機提供，采用隨體單向力模擬該載荷，作用點位于導彈質心，方向根據其推力偏斜角度設定，采用AKISPL函數擬合，得到推力曲線。

同時，武器系統發射時受到燃氣流反作用力，主要考慮燃氣流場對導流臺沖擊作用力，先通過Fluent動網格流場仿真計算得出作用力數據，然后同樣采用AKISPL函數擬合，得到發射過程燃氣沖擊作用力曲線。

2.4 整車剛柔耦合模型

整個系統模型由車體、輪胎、車軸、懸掛系統、導彈和發射臺柔性體等組成，對各結構連接依據子部件真實運動關系施加相應運動約束副。共有216個自由度。圖2為按照圖1所示拓撲連接建立的剛柔耦合動力學模型，系統的慣性坐標系如圖中所示。

3 仿真結果與分析

3.1 參數化動力學仿真

為了研究推力偏心對導彈初始擾動的影響規律，需要獲取大量的樣本數據。傳統方法上使用Adams軟件進行發射動力學仿真時，每次仿真針對的是一種固定工況，當需要進行多工況仿真時，完成一次仿真后需要人工對輸入參數進行重新設定，最后總時間花費多，同時效率極為低下。可以充分利用Adams軟件強大的二次開發功能，通過編譯能夠自動修改腳本文件(.cmd)設計變量和自動運行批處理文件(.bat)的C++程序，來實現不同設計變量條件下的自動仿真和導彈姿態樣本數據采集功能，其仿真控制流程如圖3所示。

圖3 仿真控制流程圖Fig.3 Simulation control flow chart

3.2 基于蒙特卡羅法的參數化動力學仿真結果

蒙特卡羅法[15]源于第二次世界大戰時期的德國，利用隨機數進行統計實驗，以求得的統計特征值(如均值、概率等)作為待解問題的數值解。依據相關工程經驗，推力偏斜角度隨機變量(單位(°))服從均值為0，方差為0.15的正態分布，分布區間為[-0.25，0.25]。本文采用蒙特卡羅法綜合考慮現實計算條件和時間因素，按照圖3所示的參數化自動仿真流程，以推力偏斜角度作為參數化設計變量，對其隨機采樣80次，獲得不同工況下導彈離地10 m高時姿態的多個樣本數據。

本文以決定導彈起控精度的導彈偏航方向和俯仰方向初始擾動作為研究對象，分析推力偏心對各方向初始擾動的影響。圖4所示即為對應不同推力偏斜角度情況下導彈離地10 m時的俯仰角位移、

圖4 俯仰和偏航角位移樣本數據(高度：10 m)Fig.4 Sample data of yaw and pitch angular displacement (height: 10 m)

偏航角位移樣本數據，可以看出導彈各姿態隨推力偏斜角度的不同呈現出無規律的變化。

3.3 導彈姿態分布規律分析

不同的發射裝置有著不同的力學傳遞特性，對于具有高非線性特點的發射裝置而言，分布特征明確的隨機激勵引起的初始分布規律是未知的，需要從統計學角度進行具體分析。

為了從大量樣本數據中研究對象的分布規律，通常采用直方圖法對樣本數據進行整理分析，可以繪制得到頻率直方圖。

根據導彈離地10 m時的俯仰角位移、偏航角位移樣本數據繪制所得的直方圖如圖5所示。

由各姿態的頻率直方圖形狀可以看出，在推力偏斜角度隨機變量擾動下，導彈離地10 m時俯仰角位移、偏航角位移的分布圖形近似于正態分布圖形。為了更好地了解其分布特點，假設服從正態分布，采用極大似然估計法對其概率密度函數進行參數估計。表2列出了10 m高時導彈姿態正態分布的均值和方差估計值。

圖5 俯仰和偏航角位移頻率直方圖(高度：10 m)Fig.5 Frequency histogram of pitch and yaw angular displacement (height: 10 m)

(°)

4 結束語

(1) 基于多柔體動力學理論，利用虛擬樣機技術可以快速高效地建立起復雜系統的剛柔耦合動力學模型。

(2) 本文通過編程語言對腳本文件(.cmd)的推力偏斜角度參數化設計，實現了自動仿真，得到導彈初始擾動的大量樣本數據，解決考慮隨機變量后所帶來的仿真工作量大、人工參與處理時可能導致的數據截取失誤等問題，達到了充分利用有限人力資源，提高工作效率的目的。

(3) 對于此類結構連接特性的發射系統，本文基于蒙特卡羅方法，通過參數化仿真獲取大量數據，然后運用直方圖法進行數據統計分析，研究了在推力偏斜角度隨機變量擾動下，導彈離地10 m時的俯仰角位移、偏航角位移的分布規律，得到了分布規律近似于正態分布的有益結論，并且采用極大似然估計法對其概率密度函數進行了參數估計，有助于后續起控安全性分析工作的開展。此方法對于類似隨機問題研究具有一定理論意義和參考價值。

[1] 楊錚,岳瑞華,徐中英.大型導彈發射裝置風載荷響應分析[J].現代防御技術,2015,43(5):218-222. YANG Zheng,YUE Rui-hua,XU Zhong-ying.Analysis of Wind Load Effects on Some Large Missile Launching Device[J].Modern Defence Technology,2015,43(5):218-222.

[2] 曾偉,姜毅.導彈發射車懸架半主動控制與仿真[J].現代防御技術,2015,43(1):146-150. ZENG Wei,JIANG Yi.Semi-Active Control and Simulation of Suspension System for Launch Vehicle[J].Mo-dern Defence Technology,2015,43(1):146-150.

[3] DZIOPA Z,KRZYSZTOFIK I,KORUBA Z.An Analysis of the Dynamics of A Launcher-Missile System on A Moveable Base[J].Bulletin of the Polish Academy of Sciences:Technical Sciences,2010,58(4):645-650.

[4] 陳陣,畢世華.隨機激勵對火箭初始擾動的影響[J].彈箭與制導學報,2005,25(3):70-71. CHEN Zhen,BI Shi-hua.Influence of Random Factors on Initial Disturbances[J].Journal of Projectiles,Rockets,Missiles and Guidance,2005,25(3):70-71.

[5] 李翔,畢世華,陳陣.艦船搖擺對艦載火箭初始擾動影響的多因素分析[J].北京理工大學學報,2011,31(3):253-257. LI Xiang,BI Shi-hua,CHEN Zhen.Multivariate Analysis on Initial Disturbances Effect of Ship-Launched Rocket Excited by Ship Swaying Motion[J].Transactions of Beijing Institute of Technology,2011,31(3):253-257.

[6] 康甜,賀衛亮.基于偏差模型的導彈發射動力學仿真[J].北京航空航天大學學報,2012,38(8):1111-1117. KANG Tian,HE Wei-liang.Missile Launching Dynamic Simulation Based on Variation Model[J].Journal of Beijing University of Aeronautics and Astronautics,2012,38(8):1111-1117.

[7] 張朋,譚湘霞.飛行器推力偏心及質心漂移總體估算方法研究[J].現代防御技術,2015,43(4):62-67. ZHANG Peng,TAN Xiang-xia.Estimate Methods Research on Thrust Eccentricity and Centroid Drift of the Exoatmosphere Aircraft[J].Modern Defence Technology,2015,43(4):62-67.

[8] 陳德民,槐創鋒,張克濤.精通ADAMS 2005/2007 虛擬樣機技術[M].北京:化學工業出版社,2010:174-216. CHEN De-min,HUAI Chuang-feng,ZHANG Ke-tao.Proficient in ADAMS 2005/2007 Virtual Prototyping Technology[M].Beijing:Chemical Industry Press,2010:174-216.

[9] 賀國芳.可靠性數據的收集與分析[M].北京:國防工業出版社,1995. HE Guo-fang.Reliability Data Collection and Analysis[M].Beijing:National Defense Industry Press,1995.

[10] CRAIG R R,BAMPTON M C C.Coupling of Substructures for Dynamic Analyses[J].AIAA Journal,1968,6(7):1313-1319.

[11] NEGRUT D,HARRIS B.ADAMS Theory in A Nutshell[M].US:MDI,2001:15-17.

[12] ORLANDEA N,CHANCE M A,CALAHAN D A.A Sparsity-Oriented Approach to the Dynamics Analysis and Design of Mechanical System Part[J].Journal of Engineering for Industry,1977,99(3):773-784.

[13] 洪嘉振.多體系統動力學:理論、計算方法和應用[M].上海:上海交通大學出版社,1992:28-65. HONG Jia-zhen.Dynamics of Multibody Systems:Theory,Computational Method and Application[M].Shanghai:Shanghai Jiao Tong University Press,1992:28-65.

[14] 劉德貴.動力學系統數值仿真算法[M].北京:科學出版社,2000. LIU De-gui.Digital Simulation Algorithms for Dynamic Systems[M].Beijing:Science Press,2000.

[15] DUBI A.蒙特卡洛方法在系統工程中的應用[M].衛軍胡,譯.西安:西安交通大學出版社,2007. DUBI A.Monte Carlo Applications in Systems Engineering[M].WEI Jun-hu,Translated.Xi’an:Xi’an Jiaotong University Press,2007.