趙項偉， 楊 珍， 楊 洋

(中國兵器工業(yè)試驗測試研究院，陜西 華陰 714200)

火箭橇是20世紀中后期發(fā)展起來的一種大型、高精度地面動態(tài)模擬試驗設備[1]，主要依靠靴軌配合將其約束在軌道上。隨著新式武器不斷涌現(xiàn)，火箭橇系統(tǒng)的在軌速度不斷提高[2]，系統(tǒng)運行的力學環(huán)境隨之變得苛刻，導致橇體或軌道在高速環(huán)境下面臨破壞的風險，而靴軌之間的高速沖擊是火箭橇系統(tǒng)產(chǎn)生高過載的主要原因，因此計算與分析靴軌之間的接觸特性變得尤為重要。張立乾采用ABAQUS建立橇軌火箭橇有限元模型，分析不同不平順度下的火箭橇動力特性；楊珍等[3]通過LS-DYNA計算不同工況下高超聲速火箭橇的過載均方根；Furlow[4]使用DADS建立窄軌火箭橇數(shù)值模型，分析不同參數(shù)對火箭橇響應的影響；王建[5]將火箭橇簡化為集中質量點，分析靴軌間隙、地面風速等對橇體動力學特性的影響；Lamb[6]分析了可以引發(fā)軌道共振的火箭橇臨界速度；Butova等[7]通過建立火箭橇集中質量模型計算了火箭橇系統(tǒng)在軌穩(wěn)定性。上述研究主要集中分析于不同條件下火箭橇橇體的動力學特性，對于靴軌之間接觸特性研究較少，然而在實際應用中，用于評判滑靴承載狀況的靴軌接觸力、計算靴軌接觸等效摩擦系數(shù)所需的靴軌碰撞時間、判斷軌道失效所需的靴軌接觸速度等一些關鍵參數(shù)的缺失導致火箭橇在設計過程中更多的是依靠經(jīng)驗，而非理論指導。

為了解決上述問題，本文采用Eluer-Bernouli梁單元建立橇體的有限元模型，在保留數(shù)值模型準確性的同時較全尺寸建模求解方法大幅提高求解速度，考慮軌道不平順及滑靴間隙影響獲得橇軌非線性接觸力，通過Newmark-β法結合Newton-Raphson局部迭代求解該模型，最終由靴軌接觸力、靴軌碰撞時間、靴軌接觸速度三方面入手分析靴軌之間的接觸特性，以期為火箭橇系統(tǒng)理論研究、設計及應用提供支持與指導。

1 火箭橇系統(tǒng)建模及假設

此火箭橇系統(tǒng)主要由整流罩、橇體及滑靴組成，系統(tǒng)通過固體火箭發(fā)動機推進達到設定速度，在橇體前、后滑靴及橇體質心裝有過載傳感器，在試驗中測量過載信號并通過遙測傳輸記錄。系統(tǒng)結構如圖1所示，為單軌雙滑靴結構，火箭橇系統(tǒng)在軌運行過程中，豎、側向主要受到靴軌接觸力、重力、氣動力的作用，為了簡化計算，做以下假設：

(1) 忽略滑靴磨損的影響，假設全彈道靴軌間隙不變[8-9]；

(2) 忽略橇體的滾轉效應；

(3) 軌道剛度為定值。

圖1 橇體受力示意圖Fig.1 Diagram of force on rocket sled

通過以上假設，采用Eluer-Bernouli梁單元將橇體離散為圖2所示有限元模型，靴軌接觸力、重力、氣動力施加在各節(jié)點上。

圖2 火箭橇梁單元模型Fig.2 Beam element model in rocket sled

2 含局部非線性火箭橇系統(tǒng)的動力學方程

2.1 火箭橇動力學模型

對于火箭橇系統(tǒng)，其動力學方程可寫為如下形式

(1)

針對結構外形較為簡單的火箭橇，采用梁單元便可以準確離散橇體。圖3為Eluer-Bernouli梁單元自由度示意圖[10]，每個單元由兩個節(jié)點組成，每個節(jié)點兩自由度，為x方向的位移和轉角，分別為u,θ，下標e1,e2分別表示單元的第1和第2個節(jié)點。le為單元長度，ξ為0～le中間任意一點，O為坐標原點，軸向為z向。

圖3 梁單元Fig.3 Beam element

為簡化分析，假設單元為各向同性，材料為線彈性材料且彈性模量為Ee，則：

(2)

則單元內任一截面的位移ueξ通過位移插值函數(shù)可表示為

(3)

忽略剪切變形以及旋轉慣性的影響，則Hermitian形函數(shù)Ni為

(4)

單元的應變能可寫為

(5)

將式(3)、式(4)代入式(5)可得

(6)

單元剛度矩陣Ke為

(7)

其中

(8)

單元剛度矩陣為

(9)

單元的動能為

(10)

式中，ρe為材料密度，Ae為軸段橫截面面積，且：

(11)

將式(3)、式(4)、式(11)代入式(10)可得

(12)

單元質量矩陣Me為

(13)

其中

(14)

計算得到單元質量矩陣為

(15)

加上y-z平面自由度，則單元質量矩陣和剛度矩陣可寫為

(16)

(17)

此時的qe變?yōu)?/p>

(18)

式(16)、式(17)便是梁單元的質量陣與剛度陣，橇體的阻尼則采用瑞利阻尼得到，瑞利阻尼最早由Caughey提出的廣義瑞利阻尼模型退化而來[11]，假定阻尼矩陣為剛度矩陣和質量矩陣的線性組合[12]，該理論認為總阻尼矩陣可由分塊的瑞利阻尼矩陣疊加而成，即：

C=αM+βK

(19)

α、β由下式確定

(20)

對于橇體：ξ1=ξ2=0.03，ω1與ω2為橇體的固有頻率，可通過實測獲得。

通過上述推導，便得到了橇體的M、K、C。

2.2 靴軌接觸力模型的建立

橇體的外激勵力由重力、氣動力及靴軌接觸力組成，其中重力是已知的，氣動力則可采用CFX計算若干航向速度點后插值得到，因此，重點是求出靴軌接觸力，相比于靜態(tài)接觸，橇軌力更接近于碰摩力，ADAMS內置了一種接觸算法作為描述接觸力的函數(shù)，其函數(shù)模型為[13]

(21)

其中等式右邊第一項為接觸剛度產(chǎn)生的接觸力，第二項代表接觸阻尼產(chǎn)生的阻尼力，當達到最大侵入深度后，其值不再變化，由于靴軌之間為面-面接觸，侵入深度較小，所以上式簡化為

(22)

圖4 靴軌接觸有限元模型Fig.4 Slipper-rail contact finite element model

確定接觸剛度及阻尼后，仍需確定接觸深度及速度以獲得接觸力，受限于加工裝配精度，軌道表面并非處在同一平面，而是存在高差，這種高差稱為不平順度，由于不平順度的存在，使得計算接觸深度及接觸速度時需要考慮軌道不平順度的影響，同時不平順度的取值是影響計算精度的重要因素[14-15]，為此采用實測不平順度線性插值得到計算所需不平順值。以提高數(shù)值模型準確性若滑靴位移為uei，不平順度引起的豎向位移為ε，則總侵入深度及侵入速度為

(23)

滑靴力便可寫為

(24)

考慮到靴軌之間存在間隙，以豎向為例，靴軌可有如下三種位置關系，滑靴上表面與軌道接觸、滑靴懸空、滑靴下表面與軌道接觸，因此靴軌碰撞力可寫為

(25)

式中：gap為靴軌豎向間隙;Kc1、Cc1與Kc2、Cc2分別為滑靴上、下表面的接觸剛度與接觸阻尼，側向靴軌碰撞力可同理得到。至此式(1)中M、K、C、Fc、Fg、Fa均給出表達式。

圖5 歸一化后接觸深度與速度關系Fig.5 The relation of contact depth and velocity after normalization

2.3 基于Newton-Raphson局部迭代和Newmark-β積分相結合的時域求解及驗證

采用Newmark-β法求解式(1)可得橇體響應，Newmark-β積分法是線性加速度法的一種推廣[16]，它采用如下假設

(26)

(27)

式中:Δt為時間步長;n為時間步數(shù);α和β為固定參數(shù)。一般，為了保證精度與穩(wěn)定性

(28)

上式可變形為

(29)

其中

(30)

為求解此包含非線性靴軌接觸力的動力學方程，引入Newton-Raphson局部迭代作為其非線性求解算法。

由于火箭橇側向不平順度目前并無實測數(shù)據(jù)，因此著重分析豎向橇體振動特性，由圖6各測點理論與實測過載曲線可知，橇體過載值隨航向速度增加而增加，在最大彈道速度處達到極值，理論計算包絡線與實測結果吻合較好，總體趨勢變化一致，但細觀局部可見理論與實測存在一定差異，原因可能有如下兩點：①理論計算所用軌道不平順度與實際火箭橇試驗中滑靴經(jīng)歷不平順度不一致，在每次火箭橇試驗后由于滑靴沖擊軌道不平順度均會發(fā)生改變，但由于測量軌道不平順度耗時極長，使得軌道不平順度數(shù)據(jù)更新頻率較慢；②數(shù)值模型中沒有考慮滑靴磨損所帶來的靴軌間隙變化，火箭橇在軌運行中滑靴與軌道不斷沖擊碰磨使得靴軌間隙不斷增大，進而改變靴軌接觸特性，如何準確預示滑靴磨損也是目前研究的熱點與難點。

(a) 理論彈道

(b) 質心豎向過載包絡線

(c) 前滑靴豎向過載包絡線

(d) 后滑靴豎向過載包絡線圖6 理論與實測包絡線對比Fig.6 Envelopes comparison between theory and test

3 靴軌接觸特性分析

火箭橇在軌運行中，通過滑靴將其約束在軌道中，靴軌之間的強相互作用是系統(tǒng)產(chǎn)生高過載的主因，因此，分析火箭橇在軌動力學特性的重中之重便是分析靴軌之間的耦合特性。

3.1 靴軌接觸力特性

以后滑靴為例分析靴軌接觸力特性，數(shù)值計算結果如圖7所示，靴軌接觸力隨航向速度增加而增加，最大值為2.27×106N，出現(xiàn)在2.692 s，這與航向速度最大值時刻2.687 s基本一致，同時單次碰撞分析結果表明，靴軌單次碰撞時間為0.35～0.45 ms，屬于靴軌碰撞的固有屬性，不受外因素影響，Hooser認為該時間為0.5～1 ms[17]，這可能是由于分析的橇軌類型不同導致的。

圖7 接觸力隨時間變化Fig.7 Time-domain characteristic of contact force

3.2 靴軌碰撞時間特性

火箭橇滑靴的磨損及彈道計算中，一個重要的參數(shù)便是靴軌碰撞時間，前者用其評估滑靴磨損量，后者則需要靴軌碰撞時間確定摩擦因數(shù)，為此，分析不同速度下的靴軌接觸時間占比有著重要意義。

靴軌在某一速度段下靴軌碰撞的時間與該速度段總時間的比值稱為靴軌碰撞時間占比，其隨航向速度變化曲線如圖8所示，靴軌碰撞時間占比隨航向速度增大而增加，為航向速度的正比例函數(shù)，但是由于軌道不平順隨機激勵的存在，使得靴軌碰撞時間在擬合值附近上下波動，同時前后滑靴碰撞時間占比有著巨大的差異，前滑靴在最大速度時間段靴軌碰撞時間占比達到了20.6%，而后滑靴則僅為6.8%，若定義至少有一只滑靴與軌道發(fā)生碰撞便可稱為橇軌碰撞，則由圖9可得，橇軌碰撞時間在最大速度段約占總時間的26.0%，同時由上一節(jié)靴軌單次碰撞時間基本為定值可知，隨著航向速度的增加，靴軌碰撞時間的增加是由其碰撞次數(shù)的增加而引起的。

(a) 前滑靴碰撞時間占比

圖9 火箭橇橇軌碰撞時間占比Fig.9 Sled-rail contact time proportion

3.3 靴軌接觸速度特性

火箭橇軌道的承載能力是有限的，在滑靴持續(xù)不斷沖擊下，軌道面臨斷裂的風險，現(xiàn)有研究表明，火箭橇軌道可以承受的最大沖擊速度為2.54 m/s[18]，因此，在火箭橇設計中應盡量減小靴軌沖擊速度，同時通過分析橇軌接觸速度可以確定滑軌可以承受的最大航向速度。

相比于前后滑靴的過載曲線的不同，由圖10可知前后滑靴的接觸速度波形基本一致，且接觸速度峰值均出現(xiàn)在2.692 s，這與最大接觸力出現(xiàn)時刻相同，同時前后滑靴最大接觸速度均約為2 m/s。接觸速度隨航向速度曲線變化如圖11所示，接觸速度與航向速度平方成正比，根據(jù)擬合結果，在該火箭橇設計不變情況下，若以2.54 m/s作為軌道最大承載速度，火箭橇軌道可以承受的最大速度為820 m/s，由于忽略滑靴磨損，故該值為保守值。

4 結 論

本文通過Eluer-Bernouli梁單元建立橇體有限元模型，采用Newmark-β法計算考慮靴軌耦合的火箭橇非線性模型，并通過試驗驗證了模型的正確性，分析表明：

(1) 在700 m/s速度下，橇軌接觸力峰值為2.27×106N。單次靴軌接觸時間為0.35～0.45 ms，該碰撞時間不受外界因素影響；

(2) 靴軌碰撞時間為航向速度的正比例函數(shù)，靴軌碰撞時間隨航向速度的增大而增加是由于靴軌碰撞次數(shù)增多引起的。前后滑靴碰撞時長相差巨大。在700 m/s速度下，橇軌碰撞接觸時間占總時長的26.0%；

(3) 前后滑靴接觸速度波形基本一致，接觸速度與航向速度二次方成正比例。在700 m/s航向速度下，靴軌接觸速度為2 m/s，若以2.54 m/s作為軌道可承受最大接觸速度，則該火箭橇航向速度上限保守值為820 m/s。