推力和阻力作用下柔性自旋飛行器穩定性分析

潘成龍，榮吉利，徐天富，項大林

(1.北京理工大學 宇航學院，北京 100081； 2.中國兵器工業集團航空彈藥研究院，黑龍江 哈爾濱 150030；3.北京宇航系統工程研究所，北京 100076)

0 引言

大長徑比自旋飛行器在飛行中，受到空氣動力、發動機推力、控制力等作用，會發生彈體彈性變形和振動，影響飛行器穩定性和飛行性能。高速度和高機動，使得自旋飛行器在推力作用下的穩定性問題日益突出。在研究自旋飛行器時，一般采用梁模型，很多學者對推力作用下的梁模型進行了研究。Beal[1]采用均勻梁模型，在定常推力和脈動推力作用下，分析了推力和長徑比對系統振動頻率和系統穩定性的影響。Wu[2]采用限元方法，分析了帶有集中質量的非均勻梁模型，結果表明集中質量對臨界推力值和動穩定域有相當大的影響，隨動推力作用使原剛體模態轉變成為最小的非零頻率。宋健[3]在推力和阻力作用下建立橫向運動方程，得到了均勻截面飛行器臨界推力的解析解。文獻[4-5]采用非均勻梁模型，分析了隨動推力對顫振臨界速度、穩定性及模態振型的影響。張雷等[6]考慮氣動耦合項，對推力作用下細長導彈的橫向彎曲振動進行了數學建模，討論了由氣動彈性引起的振動模態頻率以及臨界推力的變化情況。Wu等[7]采用1階活塞理論計算氣動力，分析了柔性彈箭的顫振特性。全景閣等[8]使用有限元方法分析了大長細比導彈在軸向載荷作用下的結構剛度特性。

以上研究都是基于Euler梁理論，與Euler梁相比，Timoshenko梁考慮了轉動慣量和剪切效應。陳紅永等[9]通過解析法和微分求積法(DQM)求解了梁的振動特性，分析了軸向壓力和運動效應以及軸向力導數和運動加速度對梁固有特性的影響，并對臨界載荷、臨界速度等的影響因素進行了比較研究。陳強等[10]基于動力剛度矩陣法分析了軸向變速運動彎曲梁固有頻率隨著彎曲梁軸向運動速度、加速度、軸向受力、邊界條件的變化規律。王樂等[11]采用分離變量法求解了自由邊界軸力作用下Timoshenko梁的自由振動方程。

考慮自旋時，飛行器穩定性也會受到影響。Zhu等[12]采用Timoshenko旋轉梁理論，建立彎曲、扭轉、軸向三維耦合非線性動力學方程，分析了變推力作用下系統質量偏心因素對自旋飛行器系統動力特性的影響。榮吉利等[13-15]將轉子動力學應用于外彈道學，采用有限單元法建立了隨動推力作用下柔性自旋飛行器橫向振動方程，分析了轉速和推力對系統穩定性的影響。Xu等[16]基于Timoshenko梁模型，研究了隨動推力對柔性自旋飛行器氣動彈性穩定性的影響。文獻[13-16]采用的是平均軸系，在該坐標系下系統的動量為0 kg·m/s，動量矩為0 kg·m2/s，然而其限制條件很難在真實情況滿足[17]，針對這些不足，郭東等[17]提出瞬態坐標系，陳爾康等[18]將瞬態坐標系應用于旋轉導彈動力學建模。

本文將自旋飛行器簡化為Timoshenko旋轉梁，考慮陀螺效應和剪切效應，運用轉子動力學和有限元思想，采用假設模態法和Lagrange方程，在瞬態坐標系下推導出推力和阻力作用下柔性自旋飛行器橫向振動方程。分別從轉速、推力和阻力3個方面分析了柔性自旋飛行器在軸向力作用下的穩定性，為柔性自旋飛行器研究提供了理論參考。

1 柔性彈箭動力學方程

1.1 動能、勢能和外力功

針對平均軸系法存在的不足，瞬態坐標法具有固定的坐標軸指向，能夠描述彈性變形對質心位置的影響，同時不同物理量能夠在同一坐標系中描述。

如圖1所示，建立準彈體坐標系Oxyz、瞬態坐標系Obxbybzb和微元坐標系O′ξηζ. 準彈體坐標系原點O位于未變形飛行器質心，瞬態坐標系原點Ob位于飛行器質心，兩坐標系坐標軸指向相同。忽略軸向變形，圖1中：R為軸段微元在瞬態坐標系下的位置矢量，R′為軸段微元在準彈體坐標系下的位置矢量，x0為軸段微元在準彈體坐標系下的軸向位置矢量，ΔC為飛行器質心的位置矢量，u為軸段微元中心相對準彈體坐標系的彈性變形矢量，P為隨動推力，Ω為自旋轉速。

圖1 彈體模型圖Fig.1 Flight vehicle model

在準彈體坐標系下，R′和ΔC的表達式為

R′=x0+u，

(1)

(2)

(3)

則

(4)

由(4)式知，軸段微元位置矢量在準彈體坐標系和瞬態坐標系下，相差矢量ΔC.

在準彈體坐標系下，x0和u分別為

(5)

(6)

圖2 坐標系轉換Fig.2 Transformational relation between the coordinate systems

準彈體坐標系Oxyz和微元坐標系O′ξηζ繞3個坐標系的轉換矩陣為

(7)

在微元坐標系下，軸段微元角速度[19]表示為

(8)

系統的動能可以表示為

(9)

由于彈性變形為小變形，近似取cosθη′≈1和sinθη′≈θη′≈θy，將(8)式代入(9)式，得

(10)

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

(11)

式中：EI為彎曲剛度；κGA為剪切剛度；上標符號“′”表示變量對x的偏導數。

隨動推力作的功[13]為

(12)

式中：上標符號“?”表示變量對x的3階偏導數；PN為軸向力，

(13)

隨動推力非保守部分作的虛功為

(14)

飛行器的軸向氣動力主要有頭部、尾部和底部波阻、側面摩擦阻力等，這里將軸向氣動力簡化成頭部集中阻力，集中阻力非保守部分作的虛功為[19]

δWFD=FDθy(xD,t)δuz(xD,t)-

(15)

式中：xD為頭部集中阻力作用位置。

1.2 振動方程

用有限元方法，將彈體模型沿軸向進行離散，劃分為n個梁單元。采用Timoshenko梁單元對橫向位移和轉角插值，即

(16)

式中：Φ和Ψ分別為位移和轉角插值函數矩陣；ηy和ηz分別為y軸方向和z軸方向的廣義坐標。

非保守系統Lagrange方程的形式為

(17)

(18)

式中：M、G和K分別為系統的質量矩陣、回轉矩陣和彈性剛度矩陣；KPN為軸向力保守部分得到的剛度矩陣；KP和KFD分別推力和阻力非保守部分得到的剛度矩陣；Qy和Qz分別為y軸方向和z軸方向的廣義力。矩陣和向量中的元素分別為

(19)

(20)

(21)

(22)

(23)

(24)

2 穩定性分析

將(18)式簡化為

(25)

將(25)式寫成狀態方程形式：

(26)

采用文獻[5,7]中的方法，根據李雅普諾夫穩定性理論，來判斷系統的穩定性。將系統穩定性分析轉化為求解狀態方程中矩陣Α的特征值問題。λi為非對稱矩陣Α特征值，是復數，進動頻率ωi=Im (λi)，模態阻尼為λi實部，i為階數。

對參數無量綱化[15,19]：

(27)

式中：ω為系統基頻。

以長徑比為25.3的等截面細長回轉梁模型作為仿真模型，整體分為5段，推力作用在軸段1的左端，初始時刻相關參數如表1所示，表中l為軸段長度。

表1 梁軸模型參數

2.1 平均軸系下穩定性分析

圖3 不同阻力下特征值Fig.3 Variation of eigenvalue with thrust under different drag loads

圖4所示為轉速對飛行器穩定性影響。由圖4(a)知，α=0.2下，考慮旋轉時，由于陀螺效應影響，正、反進動頻率(由于反進動頻率為負，這里取反進動絕對值進行分析)呈分離趨勢，1階反進動頻率與剛體模態頻率在P=0.8Pcr0時合并，P=1.0Pcr0時分離，1階反進動頻率增大并與2階正進動頻率合并，1階正進動頻率與2階反進動頻率在P=0.94Pcr0時合并。

圖4 不同轉速和不同阻力特征值變化Fig.4 Variation of eigenvalue with thrust under different drag loads and spinning speeds

由圖4(b)知，在(0.8Pcr0，1.0Pcr0)區間，虛部為正數，發生動態失穩，臨界推力降低，失穩區域變大，引起剛性運動和彈性振動耦合振動。

當α=1.0時，由圖4(c)和圖4(d)知，1階反進動頻率與剛體模態頻率發生耦合，2階正進動頻率逐漸減小，在(0.43Pcr0，0.76Pcr0)區間，模態頻率變為非零，靜態失穩變成動態失穩，靜態失穩推力變為動態失穩推力，失穩區域變大。

圖5所示為阻力變化對臨界推力的影響，圖中Pcr為臨界推力。由圖5可見，隨著阻力的增大，臨界推力先緩慢變小，直線下降后逐漸變小。不考慮轉速時，隨阻力增大先發生動態失穩，在α1點，動態失穩變為靜態失穩，動態臨界推力變為靜態失穩推力；考慮轉速時，隨轉速增大，臨界推力逐漸減小，當α>0.6時，轉速對臨界推力影響很小，轉速使剛體模態頻率與1階反進動頻率在α2和α3點合并，發生動態失穩。

圖5 臨界推力變化Fig.5 Changing curves of critical thrust

圖6 不同梁模型特征值變化Fig.6 Variation of eigenvalue with thrust for different beams

圖6所示為推力和阻力作用下不同梁模型穩定性的影響。將Timoshenko和Rayleigh梁模型在α=0.4時作對比分析，Timoshenko梁頻率和靜態失穩推力小于Rayleigh梁模型，對比1階頻率，剪切效應對2階頻率影響更大，剪切效應降低系統的穩定性。

當α=0和α=1.0時，計算3種均勻梁模型的臨界推力，結果如表2所示。對比表2中的數據發現，轉動慣量和剪切效應對臨界推力都有影響，與轉動慣量相比，剪切效應對臨界推力影響更大。

2.2 平均軸系和瞬態坐標系模型對比

圖7所示為轉速對飛行器穩定性影響。由圖7可見，與平均軸系模型相同，自旋轉速能夠使1階、2階正反進動頻率呈分離趨勢，對剛體模態頻率影響不大。

表2 臨界推力

圖7 瞬態坐標系中不用轉速特征值變化Fig.7 Variation of eigenvalue with thrust in the transient coordinate system at different spinning speeds

將平均軸系模型和瞬態坐標系模型進行仿真，結果如圖8～圖10所示。

圖8 不同坐標系模型特征值變化Fig.8 Variation of eigenvalue with thrust in different coordinate systems

圖9 阻力作用下不同坐標系模型特征值變化Fig.9 Variation of eigenvalue with thrust in different coordinate systems under drag load

圖10 不同模型臨界推力變化Fig.10 Variation curves of critical thrust for different models

圖8所示為轉速和阻力為0時飛行器的穩定性，可見瞬態坐標系模型產生剛體模態，發生靜態失穩，臨界推力為0.31Pcr0，與平均軸系模型相比，穩定性降低。

圖9所示為轉速0 rad/s、α=1.0時飛行器的穩定性，可見考慮阻力時，隨著推力增大，瞬態坐標系模型2階頻率不為0，剛體模態頻率和1階模態頻率合并，臨界推力為0.47Pcr0，與平均軸系模型相比，穩定性升高。

圖10所示為不同模型阻力變化對臨界推力的影響。由圖10可見，隨著阻力的增大，瞬態坐標系模型臨界推力逐漸增大，當α=0.45時，臨界推力基本不變；當α>0.81時，瞬態坐標系模型臨界推力大于平均軸系模型。

由圖8～圖10可知：阻力較小時，瞬態坐標系模型穩定性由剛體模態決定，隨著阻力增大，穩定性由剛體模態和1階彈性模態的耦合模態決定；平均軸系模型穩定性先由1階和2階彈性模態耦合模態決定，隨著阻力增大，穩定性由1階模態決定。

3 結論

本文采用Timoshenko梁模型，考慮陀螺效應和剪切效應，運用轉子動力學和有限元方法思想，建立推力作用下柔性自旋飛行器橫向振動方程，分析了阻力、轉速對柔性自旋飛行器穩定性影響。主要得到以下結論：

1)在平均軸系下：①頭部集中阻力能夠降低系統穩定性，不考慮轉速時，使1階頻率趨于0，發生靜態失穩，阻力增大，可以使2階頻率趨于0；考慮轉速時，自旋能夠使1階反進動頻率與剛體模態頻率在合并，發生動態失穩，失穩區域變大。②阻力能夠使臨界推力和臨界轉速降低，轉速使靜態失穩變為動態失穩。③轉動慣量和剪切效應對穩定性都有影響，與轉動慣量相比，剪切效應影響更大，特別是對2階頻率影響。

2)在瞬態坐標系下：推力作用產生剛體模態，系統穩定性主要由剛體模態決定，阻力作用能提高系統穩定性；自旋轉速不改變失穩區域。