車載導(dǎo)彈推力偏心初始擾動影響研究*

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

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

0 引言

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

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

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

1.1 修正的Craig-Bampton

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

易知，在C-B模態(tài)坐標(biāo)下，子結(jié)構(gòu)的結(jié)構(gòu)動力學(xué)方程為

(1)

式中：Mi,Ci,Ki,Ri分別為C-B模態(tài)基下的子結(jié)構(gòu)質(zhì)量、阻尼、剛度和外力矩陣。

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

(2)

(3)

因此，子結(jié)構(gòu)的物理自由度可以表示為

(4)

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

1.2 多體系統(tǒng)剛?cè)狁詈蟿恿W(xué)方程

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

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

(5)

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

將柔性體的計算結(jié)果和多剛體理論方法相結(jié)合可以進(jìn)一步得到多體系統(tǒng)的剛?cè)狁詈蟿恿W(xué)方程。通過文獻(xiàn)[12-14]可以推導(dǎo)出使用拉格朗日乘子得到的第i個柔體或剛體的方程形式為

(6)

式中：K為動能的表達(dá)式;Qi為廣義力，包括單元彈性變形和外加載荷(包括單元間載荷)引起的廣義力。將式(11)與約束方程C(q,t)=0聯(lián)立，即構(gòu)成剛?cè)狁詈隙囿w動力學(xué)方程。

2 車載導(dǎo)彈剛?cè)狁詈蟿恿W(xué)模型

2.1 武器系統(tǒng)組成分析與假設(shè)

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

根據(jù)車載導(dǎo)彈武器系統(tǒng)的子部件功能關(guān)系，對整車系統(tǒng)進(jìn)行合理處理，簡化后的模型主要分為車輪、調(diào)平支腿、懸掛系統(tǒng)、車體、上裝設(shè)備、導(dǎo)流臺、發(fā)射臺、起豎油缸和導(dǎo)彈。其中，導(dǎo)彈無約束豎立在發(fā)射臺支撐盤上；發(fā)射臺通過旋轉(zhuǎn)鉸與車體和起豎油缸連接；導(dǎo)流臺、車輪和調(diào)平支腿均與地面接觸；上裝設(shè)備緊固在車體上；車體與車輪通過懸掛系統(tǒng)連接。系統(tǒng)主要部分具體拓?fù)溥B接如圖1所示。

為了描述武器系統(tǒng)垂直發(fā)射時的動力學(xué)特性，特作如下假設(shè)：

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

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

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

(4) 不考慮發(fā)射過程中導(dǎo)彈的變質(zhì)量特性。

2.2 柔性體建模

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

圖1 拓?fù)溥B接圖Fig.1 Topological connection graph

模態(tài)階數(shù)789101112模態(tài)頻率/Hz35.943.252.562.880.7102.8

2.3 激勵

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

導(dǎo)彈推力由發(fā)動機提供，采用隨體單向力模擬該載荷，作用點位于導(dǎo)彈質(zhì)心，方向根據(jù)其推力偏斜角度設(shè)定，采用AKISPL函數(shù)擬合，得到推力曲線。

同時，武器系統(tǒng)發(fā)射時受到燃?xì)饬鞣醋饔昧?，主要考慮燃?xì)饬鲌鰧?dǎo)流臺沖擊作用力，先通過Fluent動網(wǎng)格流場仿真計算得出作用力數(shù)據(jù)，然后同樣采用AKISPL函數(shù)擬合，得到發(fā)射過程燃?xì)鉀_擊作用力曲線。

2.4 整車剛?cè)狁詈夏Ｐ?/h3>

整個系統(tǒng)模型由車體、輪胎、車軸、懸掛系統(tǒng)、導(dǎo)彈和發(fā)射臺柔性體等組成，對各結(jié)構(gòu)連接依據(jù)子部件真實運動關(guān)系施加相應(yīng)運動約束副。共有216個自由度。圖2為按照圖1所示拓?fù)溥B接建立的剛?cè)狁詈蟿恿W(xué)模型，系統(tǒng)的慣性坐標(biāo)系如圖中所示。

3 仿真結(jié)果與分析

3.1 參數(shù)化動力學(xué)仿真

為了研究推力偏心對導(dǎo)彈初始擾動的影響規(guī)律，需要獲取大量的樣本數(shù)據(jù)。傳統(tǒng)方法上使用Adams軟件進(jìn)行發(fā)射動力學(xué)仿真時，每次仿真針對的是一種固定工況，當(dāng)需要進(jìn)行多工況仿真時，完成一次仿真后需要人工對輸入?yún)?shù)進(jìn)行重新設(shè)定，最后總時間花費多，同時效率極為低下?？梢猿浞掷肁dams軟件強大的二次開發(fā)功能，通過編譯能夠自動修改腳本文件(.cmd)設(shè)計變量和自動運行批處理文件(.bat)的C++程序，來實現(xiàn)不同設(shè)計變量條件下的自動仿真和導(dǎo)彈姿態(tài)樣本數(shù)據(jù)采集功能，其仿真控制流程如圖3所示。

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

3.2 基于蒙特卡羅法的參數(shù)化動力學(xué)仿真結(jié)果

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

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

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

偏航角位移樣本數(shù)據(jù)，可以看出導(dǎo)彈各姿態(tài)隨推力偏斜角度的不同呈現(xiàn)出無規(guī)律的變化。

3.3 導(dǎo)彈姿態(tài)分布規(guī)律分析

不同的發(fā)射裝置有著不同的力學(xué)傳遞特性，對于具有高非線性特點的發(fā)射裝置而言，分布特征明確的隨機激勵引起的初始分布規(guī)律是未知的，需要從統(tǒng)計學(xué)角度進(jìn)行具體分析。

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

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

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

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

(°)

4 結(jié)束語

(1) 基于多柔體動力學(xué)理論，利用虛擬樣機技術(shù)可以快速高效地建立起復(fù)雜系統(tǒng)的剛?cè)狁詈蟿恿W(xué)模型。

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

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

[1] 楊錚,岳瑞華,徐中英.大型導(dǎo)彈發(fā)射裝置風(fēng)載荷響應(yīng)分析[J].現(xiàn)代防御技術(shù),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] 曾偉,姜毅.導(dǎo)彈發(fā)射車懸架半主動控制與仿真[J].現(xiàn)代防御技術(shù),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].彈箭與制導(dǎo)學(xué)報,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].北京理工大學(xué)學(xué)報,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] 康甜,賀衛(wèi)亮.基于偏差模型的導(dǎo)彈發(fā)射動力學(xué)仿真[J].北京航空航天大學(xué)學(xué)報,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] 張朋,譚湘霞.飛行器推力偏心及質(zhì)心漂移總體估算方法研究[J].現(xiàn)代防御技術(shù),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] 陳德民,槐創(chuàng)鋒,張克濤.精通ADAMS 2005/2007 虛擬樣機技術(shù)[M].北京:化學(xué)工業(yè)出版社,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] 賀國芳.可靠性數(shù)據(jù)的收集與分析[M].北京:國防工業(yè)出版社,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] 洪嘉振.多體系統(tǒng)動力學(xué):理論、計算方法和應(yīng)用[M].上海:上海交通大學(xué)出版社,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] 劉德貴.動力學(xué)系統(tǒng)數(shù)值仿真算法[M].北京:科學(xué)出版社,2000. LIU De-gui.Digital Simulation Algorithms for Dynamic Systems[M].Beijing:Science Press,2000.

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