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

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

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

0 引言

現代火箭或導彈等空間飛行器通常選取較高的推重比，結構長徑比設計也越來越大，其目的是為了提高飛行速度及增加有效射程。但這同時也帶來不利影響，導致彈箭結構的彈性效應愈加明顯，柔性飛行器的結構振動特性和動力學穩定性問題逐漸得到學者們的關注。

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

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

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

1 動能、勢能和外力功

1.1 動能和勢能

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

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

圖2 坐標系轉換關系Fig.2 The transformational relation between the coordinate systems

忽略彈體的軸向變形，即假設結構無縱向伸縮，但是考慮彎曲和剪切作用對軸向變形的影響。系統動能為

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

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

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

考慮剪切變形的影響，系統的彈性勢能為

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

1.2 外力功

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

同樣地，得到Oyz 平面內的軸向位移

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

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

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

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

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

2 運動方程

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

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

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

式中:

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

非保守系統的Lagrange 方程表達為

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

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

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

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

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

3 數值仿真及結果分析

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

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

表1 梁軸模型參數Tab.1 Parameters of beam model

表2 臨界轉速Tab.2 Critical speed

表3 臨界推力Tab.3 Critical thrust

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

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

圖4 不同條件下系統質心的位移響應Fig.4 Displacement responses of mass center under different conditions

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

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

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

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

圖8 不同轉速下系統質心的位移響應Fig.8 Displacement responses of mass center under different spinning speeds

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

圖9 臨界轉速變化曲線Fig.9 Variation curves of critical spin speed

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

4 結論

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

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

2)推力達到或超過系統臨界推力時，系統發生失穩;轉速的增加會降低系統的臨界推力。

3)模型的非均勻性使結構的慣性、剛度和軸向力分布等發生變化，從而對系統的臨界轉速和臨界推力造成很大影響。

綜上，彈體自身旋轉和隨動推力的共同作用使系統的振動特性發生改變，在飛行器設計時需要引起關注。

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］許赟，謝長川，楊超. 推力作用下細長彈箭橫向振動及穩定性分析［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 梁模型的旋轉彈箭橫向振動模態分析［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］榮吉利，李瑞英. 大長徑比自旋彈箭橫向自振特性的有限元計算方法與結果分析［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］榮吉利，徐天富，王璽，等. 隨動推力作用下柔性旋轉飛行器動力學穩定性分析［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)