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

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

0 引言

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

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

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

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

1 動(dòng)能、勢(shì)能和外力功

1.1 動(dòng)能和勢(shì)能

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

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

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

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

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

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

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

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

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

1.2 外力功

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

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

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

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

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

所以隨動(dòng)推力保守部分所做功為

在彈體尾部，隨動(dòng)推力非保守部分所做虛功為

2 運(yùn)動(dòng)方程

按有限元方法的思想，將彈體模型沿軸向進(jìn)行離散，整體劃分為n 個(gè)梁?jiǎn)卧?。為了便于分析，每個(gè)梁?jiǎn)卧茷榈冉孛媪?。將整體坐標(biāo)系轉(zhuǎn)換到單元局部坐標(biāo)系下，采用Timoshenko 梁?jiǎn)卧獙?duì)其橫向位移和轉(zhuǎn)角進(jìn)行插值，得

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

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

式中:

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

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

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

式中:Ms、ΩJs和Ks分別為單元的質(zhì)量矩陣、回轉(zhuǎn)矩陣和剛度矩陣，此時(shí)廣義力Q1s和Q2s為不包括隨動(dòng)推力在內(nèi)的廣義力。單元?jiǎng)偠染仃嚢穗S動(dòng)推力的影響，Ks=Kes-Kcs-Kncs，Kes為單元彈性剛度矩陣，Kcs為由推力保守部分得到的單元?jiǎng)偠染仃?，Kncs為由推力非保守部分得到的單元?jiǎng)偠染仃?。單元矩陣Ms、Js和Kes的具體表述請(qǐng)參考文獻(xiàn)［10］，這里只給出單元矩陣Kcs和Kncs的表述:

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

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

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

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

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

以長(zhǎng)徑比為25 的等截面細(xì)長(zhǎng)回轉(zhuǎn)梁作為彈體仿真模型，整體分為4 段，相關(guān)參數(shù)見(jiàn)表1. 在不同條件下對(duì)此模型分別進(jìn)行臨界轉(zhuǎn)速［10-11］和臨界推力分析［12］，結(jié)果分別見(jiàn)表2和表3. 表2中:無(wú)量綱推力=P/Pcr0，Pcr0為Ω=0 時(shí)系統(tǒng)的臨界推力，Pcr0=9.84 ×106N;無(wú)量綱臨界轉(zhuǎn)速= ωcr/ωcr0. 表3中:無(wú)量綱臨界推力=Pcr/Pcr0;無(wú)量綱轉(zhuǎn)速=Ω/ωcr0，ωcr0為P =0 時(shí)系統(tǒng)的一階臨界轉(zhuǎn)速，ωcr0=97.16 rad/s. 可見(jiàn)，不同條件下系統(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)橫向振動(dòng)響應(yīng)問(wèn)題，把偏心力作為彈體結(jié)構(gòu)的激振力施加于彈體質(zhì)心位置處，設(shè)偏心距為e =1 mm，激振力頻率與自轉(zhuǎn)頻率相同，即計(jì)算轉(zhuǎn)速分別為0.5、1.0、1.5 時(shí)和推力分別為0、0.1 時(shí)系統(tǒng)質(zhì)心的橫向位移響應(yīng)情況，結(jié)果如圖4所示。

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

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

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

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

圖6 臨界推力作用下系統(tǒng)質(zhì)心的位移響應(yīng)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ì)心的位移響應(yīng)Fig.8 Displacement responses of mass center under different spinning speeds

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

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

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

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

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

綜上，彈體自身旋轉(zhuǎn)和隨動(dòng)推力的共同作用使系統(tǒng)的振動(dòng)特性發(fā)生改變，在飛行器設(shè)計(jì)時(shí)需要引起關(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］王良明. 大長(zhǎng)徑比彈箭在飛行時(shí)的柔性變形特性分析［J］. 兵工學(xué)報(bào)，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］王良明，王中原，易文俊. 大長(zhǎng)徑比彈箭飛行中的共振條件研究［J］. 彈道學(xué)報(bào)，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］許赟，謝長(zhǎng)川，楊超. 推力作用下細(xì)長(zhǎng)彈箭橫向振動(dòng)及穩(wěn)定性分析［J］. 工程力學(xué)，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］張雷，王永. 尾部有推力自由梁的振動(dòng)分析［J］. 彈道學(xué)報(bào)，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］何斌，芮筱亭，陸毓琪. 柔性彈箭飛行力學(xué)建模研究［J］. 彈道學(xué)報(bào)，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)彈箭橫向振動(dòng)模態(tài)分析［J］. 北京理工大學(xué)學(xué)報(bào)，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］榮吉利，李瑞英. 大長(zhǎng)徑比自旋彈箭橫向自振特性的有限元計(jì)算方法與結(jié)果分析［J］. 兵工學(xué)報(bào)，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］榮吉利，徐天富，王璽，等. 隨動(dòng)推力作用下柔性旋轉(zhuǎn)飛行器動(dòng)力學(xué)穩(wěn)定性分析［J］. 宇航學(xué)報(bào)，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)