隨動推力作用下柔性自旋飛行器橫向振動響應及失穩(wěn)分析

榮吉利，徐天富，王璽，殷新喆

(1.北京理工大學宇航學院，北京100081;2. 北京航天發(fā)射技術(shù)研究所，北京100076)

0 引言

現(xiàn)代火箭或?qū)椀瓤臻g飛行器通常選取較高的推重比，結(jié)構(gòu)長徑比設計也越來越大，其目的是為了提高飛行速度及增加有效射程。但這同時也帶來不利影響，導致彈箭結(jié)構(gòu)的彈性效應愈加明顯，柔性飛行器的結(jié)構(gòu)振動特性和動力學穩(wěn)定性問題逐漸得到學者們的關(guān)注。

Chae 等［1］考慮了多種因素對細長柔性飛行器進行了動力分析以及氣動分析，得到了臨界推力作用下的發(fā)散響應結(jié)果。Trikha 等［2］采用兩種方法推導了柔性空間飛行器的一般運動方程，其分析結(jié)果表明速度、加速度、隨動推力等因素對結(jié)構(gòu)穩(wěn)定性產(chǎn)生重要影響。文獻［3 -4］提出了大長徑比彈箭柔性彎曲的動力學模型和有限元模型，并研究了彈箭飛行過程中的共振不穩(wěn)定性問題，得出柔性變形使彈箭的飛行穩(wěn)定性變壞的結(jié)論。許赟等［5］建立了隨動推力作用下的非均勻梁振動分析數(shù)學模型，通過數(shù)值仿真表明隨動推力會誘發(fā)彈箭飛行器的動力學失穩(wěn)，同時影響結(jié)構(gòu)的振動頻率和振型特性。張雷等［6］引入氣動耦合項對有推力的細長導彈的橫向彎曲振動進行了數(shù)學建模，其分析結(jié)果表明氣動彈性使得結(jié)構(gòu)穩(wěn)定性降低。以上這些研究都是針對非自旋類柔性飛行器的，都沒有考慮飛行器自旋的作用。

何斌等［7］采用柔體浮動框架和氣動彈性小參數(shù)攝動法建立了柔性旋轉(zhuǎn)彈箭飛行力學模型，給出了被動段時系統(tǒng)臨界狀態(tài)的條件和臨界轉(zhuǎn)速。Zhu等［8］考察了變推力作用下系統(tǒng)質(zhì)量偏心因素對自旋飛行器系統(tǒng)動力特性的影響，其仿真結(jié)果表明結(jié)構(gòu)振動規(guī)律與變推力頻率密切相關(guān)。目前以柔性自旋飛行器為對象考慮隨動推力作用的研究還較少。本文在前人工作的基礎上，著重考慮飛行器的自身旋轉(zhuǎn)和隨動推力的共同作用，對飛行器橫向振動響應及失穩(wěn)情況進行分析。

本文采用剪切變形對軸向位移有貢獻的Timoshenko 梁模型，基于有限元方法建立計入陀螺力矩及隨動推力影響的系統(tǒng)運動方程，通過數(shù)值仿真分析了一定轉(zhuǎn)速時隨動推力作用下系統(tǒng)的橫向振動響應及結(jié)構(gòu)失穩(wěn)情況。

1 動能、勢能和外力功

1.1 動能和勢能

柔性自旋飛行器的簡化模型如圖1所示，圖中Oxyz 為準彈體坐標系，彈體以角速度Ω 繞x 軸自旋;隨體坐標系O'ξηζ 固連在軸段微元形心上，跟隨彈體微元一起轉(zhuǎn)動。如圖2所示，由兩坐標系之間的轉(zhuǎn)換關(guān)系，軸段微元角速度在隨體坐標系下的分量可表示為

圖1 彈體模型Fig.1 Model of flight vehicle

圖2 坐標系轉(zhuǎn)換關(guān)系Fig.2 The transformational relation between the coordinate systems

忽略彈體的軸向變形，即假設結(jié)構(gòu)無縱向伸縮，但是考慮彎曲和剪切作用對軸向變形的影響。系統(tǒng)動能為

式中:uy、uz分別為梁截面y 向、z 向的橫向位移;μ、jd和jp分別為單位長度的質(zhì)量密度、直徑轉(zhuǎn)動慣量和極轉(zhuǎn)動慣量;l 為彈體長度。

當彈體彈性變形較小時，可取近似關(guān)系cos θη'≈1和sin θη'≈θη'≈θy，略去高階小量，得到簡化動能表達式

式中:θy、θz分別為該截面的轉(zhuǎn)角;-2jpΩθ·zθy為耦合項，是由陀螺力矩作用引起的。

考慮剪切變形的影響，系統(tǒng)的彈性勢能為

式中:EI 和κGA 分別為彎曲剛度和剪切剛度。

1.2 外力功

考慮剪切變形后，其將與橫向彎曲位移共同對軸向位移產(chǎn)生影響［9］，圖3所示為Oxy 平面內(nèi)在剪切力Fsy作用下微元的軸向位移示意圖，由圖中關(guān)系可得Oxy 平面內(nèi)的軸向位移

同樣地，得到Oyz 平面內(nèi)的軸向位移

式中:γy為微元中心軸線繞y 軸的剪切角。

忽略氣動力作用，將推力P 視為定常隨動推力，在推力和慣性力作用下得到軸向力分布表達式

圖3 Timoshenko 梁的軸向位移Fig.3 Axial displacement of Timoshenko beam

所以隨動推力保守部分所做功為

在彈體尾部，隨動推力非保守部分所做虛功為

2 運動方程

按有限元方法的思想，將彈體模型沿軸向進行離散，整體劃分為n 個梁單元。為了便于分析，每個梁單元近似為等截面梁。將整體坐標系轉(zhuǎn)換到單元局部坐標系下，采用Timoshenko 梁單元對其橫向位移和轉(zhuǎn)角進行插值，得

式中:N 和D 分別為位移和轉(zhuǎn)角插值函數(shù)矩陣;u1s和u2s分別為y 向和z 向的廣義坐標［10］。

離散化后第i 個單元的軸向力可表示為

式中:

其中msj為第j 個單元的質(zhì)量，j =1，2，…，i -1，μi為第i 個單元的線密度，s 為單元局部坐標。

非保守系統(tǒng)的Lagrange 方程表達為

式中:q 和Q 分別為廣義坐標和包括非保守力在內(nèi)的廣義力。取q 分別為u1s和u2s，由此得到單元的運動方程

式中:Ms、ΩJs和Ks分別為單元的質(zhì)量矩陣、回轉(zhuǎn)矩陣和剛度矩陣，此時廣義力Q1s和Q2s為不包括隨動推力在內(nèi)的廣義力。單元剛度矩陣包含了隨動推力的影響，Ks=Kes-Kcs-Kncs，Kes為單元彈性剛度矩陣，Kcs為由推力保守部分得到的單元剛度矩陣，Kncs為由推力非保守部分得到的單元剛度矩陣。單元矩陣Ms、Js和Kes的具體表述請參考文獻［10］，這里只給出單元矩陣Kcs和Kncs的表述:

式中:li為單元長度。注意，(16)式表述的Kncs只作用在隨動推力所施加的單元上，在其他單元的剛度矩陣中Kncs元素皆為0.

考慮阻尼影響，利用有限元方法，可以得到柔性自旋飛行器運動方程的一般形式，即

式中:M、C、G 和K 分別為系統(tǒng)的質(zhì)量矩陣、阻尼矩陣、回轉(zhuǎn)矩陣和剛度矩陣;Q 為作用在彈體上載荷列陣。陀螺力矩的作用使得回轉(zhuǎn)矩陣G 為反對稱矩陣，而隨動推力的作用使得剛度矩陣K 為非對稱矩陣。最后施加平均彈軸條件［3］，對自由邊界的彈體變形進行約束。

3 數(shù)值仿真及結(jié)果分析

根據(jù)前文得到的細長自旋飛行器系統(tǒng)的運動方程，采用Newmark 方法進行求解，編制了相應的仿真程序，對在推力和自旋作用下的彈體振動情況進行仿真分析。

以長徑比為25 的等截面細長回轉(zhuǎn)梁作為彈體仿真模型，整體分為4 段，相關(guān)參數(shù)見表1. 在不同條件下對此模型分別進行臨界轉(zhuǎn)速［10-11］和臨界推力分析［12］，結(jié)果分別見表2和表3. 表2中:無量綱推力=P/Pcr0，Pcr0為Ω=0 時系統(tǒng)的臨界推力，Pcr0=9.84 ×106N;無量綱臨界轉(zhuǎn)速= ωcr/ωcr0. 表3中:無量綱臨界推力=Pcr/Pcr0;無量綱轉(zhuǎn)速=Ω/ωcr0，ωcr0為P =0 時系統(tǒng)的一階臨界轉(zhuǎn)速，ωcr0=97.16 rad/s. 可見，不同條件下系統(tǒng)的臨界轉(zhuǎn)速和臨界推力是不同的。

表1 梁軸模型參數(shù)Tab.1 Parameters of beam model

表2 臨界轉(zhuǎn)速Tab.2 Critical speed

表3 臨界推力Tab.3 Critical thrust

分析由質(zhì)量偏心因素引起的彈箭結(jié)構(gòu)橫向振動響應問題，把偏心力作為彈體結(jié)構(gòu)的激振力施加于彈體質(zhì)心位置處，設偏心距為e =1 mm，激振力頻率與自轉(zhuǎn)頻率相同，即計算轉(zhuǎn)速分別為0.5、1.0、1.5 時和推力分別為0、0.1 時系統(tǒng)質(zhì)心的橫向位移響應情況，結(jié)果如圖4所示。

由于隨動推力的作用，使得彈體結(jié)構(gòu)的彎曲和剪切變形對軸向位移產(chǎn)生影響，進而減小了系統(tǒng)剛度，所以推力的增加會使結(jié)構(gòu)的振動響應幅值增大。由圖4(a)可見，仿真結(jié)果與理論相符。圖4(b)中結(jié)果顯示，曲線1 和曲線3 兩種情況表明系統(tǒng)發(fā)生共振，結(jié)構(gòu)變形發(fā)散，因為它們的轉(zhuǎn)速均達到了各自狀況下的臨界轉(zhuǎn)速(見表2)。圖4(b)中曲線2 情況表明結(jié)構(gòu)產(chǎn)生拍振現(xiàn)象，原因是隨動推力的增加降低了系統(tǒng)的臨界轉(zhuǎn)速，使此時的轉(zhuǎn)速避開了本狀況下的臨界轉(zhuǎn)速，因此也避免了共振的發(fā)生。對應于圖4(b)中曲線1、曲線3 兩種情況，圖5顯示了彈體質(zhì)心橫向位移在臨界轉(zhuǎn)速下的共振發(fā)散軌跡。對比圖4(a)和圖4(c)，在相同推力作用下，轉(zhuǎn)速的增加使質(zhì)量偏心力增大，從而使得相應的振動響應幅值增大。

圖4 不同條件下系統(tǒng)質(zhì)心的位移響應Fig.4 Displacement responses of mass center under different conditions

圖5 彈體質(zhì)心軌跡Fig.5 Trace of mass center of flight vehicle

現(xiàn)考察推力接近或超過臨界推力時彈體結(jié)構(gòu)的響應問題。圖6顯示了P=1.00 時在不同轉(zhuǎn)速下系統(tǒng)的位移響應，明顯發(fā)現(xiàn)此時不論轉(zhuǎn)速如何系統(tǒng)都將產(chǎn)生失穩(wěn)。因為彈體自旋在一定程度上降低了系統(tǒng)的臨界推力，使得此時的推力Pcr0大于任何轉(zhuǎn)速下的臨界推力(見表3)，所以不論轉(zhuǎn)速如何最終都將導致系統(tǒng)失穩(wěn)。圖7顯示了彈體顫振失穩(wěn)時質(zhì)心的發(fā)散軌跡。事實上，當推力達到臨界推力值后，由于結(jié)構(gòu)振動失穩(wěn)，圖6和圖7中的情況會造成系統(tǒng)結(jié)構(gòu)的破壞，而實際位移不會如此之大。取推力=0.76，不同轉(zhuǎn)速下的系統(tǒng)質(zhì)心位移響應如圖8所示。對于圖8中曲線1、曲線2 兩種情況，推力均未達到各自狀況下的臨界推力(見表3)，所以振動響應曲線穩(wěn)定;對于曲線3 情況，由于彈體自旋作用使得此時的臨界推力值降幅很大(見表3)，推力已達到臨界推力值，所以結(jié)構(gòu)發(fā)生顫振失穩(wěn)。

圖6 臨界推力作用下系統(tǒng)質(zhì)心的位移響應Fig.6 Displacement response of mass center under critical follower thrust

圖7 彈體質(zhì)心軌跡Fig.7 Trace of mass center of flight vehicle

圖8 不同轉(zhuǎn)速下系統(tǒng)質(zhì)心的位移響應Fig.8 Displacement responses of mass center under different spinning speeds

圖9為不同梁軸模型的前兩階臨界轉(zhuǎn)速隨推力增加而變化的曲線，其變化規(guī)律與推力作用下的進動頻率變化規(guī)律相似，發(fā)生了頻率合并的現(xiàn)象。可以發(fā)現(xiàn)，隨著推力的增加，不同模型的一階無量綱臨界轉(zhuǎn)速的差異也逐漸增大。圖10 為不同梁軸模型的臨界推力隨轉(zhuǎn)速增加而變化的曲線。可見，均勻梁軸模型的陀螺效應很小，對臨界推力的影響可以忽略;而對于非均勻梁軸模型，軸向力的分布發(fā)生改變，同時陀螺效應隨轉(zhuǎn)速的增加而增大，當轉(zhuǎn)速達到某一值時，使得系統(tǒng)的剛性模態(tài)頻率與彈性模態(tài)頻率產(chǎn)生合并，從而使系統(tǒng)的臨界推力突然減小。此情況對應于圖10 中的階躍部分，同時與表3數(shù)據(jù)對應，圖8中的振動情況也能得以解釋。所以，結(jié)構(gòu)模型的非均勻性對系統(tǒng)的臨界轉(zhuǎn)速以及臨界推力的影響都非常大。

圖9 臨界轉(zhuǎn)速變化曲線Fig.9 Variation curves of critical spin speed

圖10 臨界推力變化曲線Fig.10 Variation curves of critical thrust

4 結(jié)論

以柔性自旋飛行器為研究對象，考慮了剪切變形對軸向位移的影響，采用有限元方法對系統(tǒng)進行了離散，推導出含有陀螺力矩及隨動推力影響的系統(tǒng)運動方程，通過算例仿真對轉(zhuǎn)速和隨動力推力作用下的系統(tǒng)橫向振動響應和結(jié)構(gòu)失穩(wěn)情況進行了分析，得到以下結(jié)論:

1)隨動推力的增加會減小系統(tǒng)剛度，使振動幅值增大，而且會減小自旋飛行器的臨界轉(zhuǎn)速;當激振頻率等于系統(tǒng)臨界轉(zhuǎn)速頻率時，系統(tǒng)產(chǎn)生共振。

2)推力達到或超過系統(tǒng)臨界推力時，系統(tǒng)發(fā)生失穩(wěn);轉(zhuǎn)速的增加會降低系統(tǒng)的臨界推力。

3)模型的非均勻性使結(jié)構(gòu)的慣性、剛度和軸向力分布等發(fā)生變化，從而對系統(tǒng)的臨界轉(zhuǎn)速和臨界推力造成很大影響。

綜上，彈體自身旋轉(zhuǎn)和隨動推力的共同作用使系統(tǒng)的振動特性發(fā)生改變，在飛行器設計時需要引起關(guān)注。

References)

［1］Chae S，Hodges D H. Dynamics and aeroelastic analysis of missiles［C］∥Proceedings of the 44th AIAA/ASME/ASCE/AHS Structures，Structural Dynamics，and Materials Conference. Norfolk，VA:SIAM，2003:1 -6.

［2］Trikha M，Mahapatra D R，Gopalakrishnan S，et al. Structural stability of slender aerospace vehicles:part Ⅱ:numerical simulations［J］. International Journal of Mechanical Sciences，2010，52(9):1145 -1157.

［3］王良明. 大長徑比彈箭在飛行時的柔性變形特性分析［J］. 兵工學報，2000，21(2):108 -111.WANG Liang-ming. An analysis on the flexibility in flight of projectiles or rockets having high L/D ratios［J］. Acta Armamentarii，2000，21(2):108 -111.(in Chinese)

［4］王良明，王中原，易文俊. 大長徑比彈箭飛行中的共振條件研究［J］. 彈道學報，2001，13(1):51 -54.WANG Liang-ming，WANG Zhong-yuan，YI Wen-jun. Resonance unstability in ballistics of flexible body［J］. Journal of Ballistics，2001，13(1):51 -54.(in Chinese)

［5］許赟，謝長川，楊超. 推力作用下細長彈箭橫向振動及穩(wěn)定性分析［J］. 工程力學，2009，26(12):211 -215.XU Yun，XIE Chang-chuan，YANG Chao. Transverse vibration and dynamic stability analysis of slender projects under thrust［J］.Engineering Mechanics，2009，26(12):211 -215.(in Chinese)

［6］張雷，王永. 尾部有推力自由梁的振動分析［J］. 彈道學報，2010，22(1):83 -86.ZHANG Lei，WANG Yong. Vibration analysis of free-free beam with end rocket thrust［J］. Journal of Ballistics，2010，22(1):83-86.(in Chinese)

［7］何斌，芮筱亭，陸毓琪. 柔性彈箭飛行力學建模研究［J］. 彈道學報，2006，18(1):22 -24.HE Bin，RUI Xiao-ting，LU Yu-qi. A study on flight dynamic modeling of flexible shell/rocket［J］. Journal of Ballistics，2006，18(1):22 -24.(in Chinese)

［8］Zhu H L，Yu L，Liang S H，et al. Dynamic simulation of a flexible spinning vehicle with the periodic varing thrust［J］. Journal of Shanghai University:English Edition，2008，12(2):115 -119.

［9］Choi S H，Pierre C，Ulsoy A G. Consistent modeling of rotating timoshenko shafts subject to axial loads［J］. Journal of Vibration and Acoustics-Transactions of the Asme，1992，114(2):249 -259.

［10］榮吉利，徐天富，李健. 基于Timoshenko 梁模型的旋轉(zhuǎn)彈箭橫向振動模態(tài)分析［J］. 北京理工大學學報，2012，32(6):551 -555.RONG Ji-li，XU Tian-fu，LI Jian. Modal analysis of transverse vibration of a spinning rocket based on Timoshenko beam［J］.Transactions of Beijing Institute of Technology，2012，32(6):551 -555.(in Chinese)

［11］榮吉利，李瑞英. 大長徑比自旋彈箭橫向自振特性的有限元計算方法與結(jié)果分析［J］. 兵工學報，2002，23(1):79 -82.RONG Ji-li，LI Rui-ying. Finite element computational method and resultant analysis of the transverse self-oscillation characteristics of a slender spinning rocket［J］. Acta Armamentarii，2002，23(1):79 -82.(in Chinese)

［12］榮吉利，徐天富，王璽，等. 隨動推力作用下柔性旋轉(zhuǎn)飛行器動力學穩(wěn)定性分析［J］. 宇航學報，2015，36(1):18 -24.RONG Ji-li，XU Tian-fu，WANG Xi，et al. Dynamic stability analysis of slender spinning flight vehicles under follower thrust［J］. Journal of Astronautics，2015，36(1):18 -24. (in Chinese)