推力和阻力作用下柔性自旋飛行器穩(wěn)定性分析

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

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

0 引言

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

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

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

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

1 柔性彈箭動力學方程

1.1 動能、勢能和外力功

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

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

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

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

R′=x0+u，

(1)

(2)

(3)

則

(4)

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

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

(5)

(6)

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

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

(7)

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

(8)

系統(tǒng)的動能可以表示為

(9)

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

(10)

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

(11)

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

隨動推力作的功[13]為

(12)

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

(13)

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

(14)

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

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

(15)

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

1.2 振動方程

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

(16)

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

非保守系統(tǒng)Lagrange方程的形式為

(17)

(18)

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

(19)

(20)

(21)

(22)

(23)

(24)

2 穩(wěn)定性分析

將(18)式簡化為

(25)

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

(26)

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

對參數(shù)無量綱化[15,19]：

(27)

式中：ω為系統(tǒng)基頻。

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

表1 梁軸模型參數(shù)

2.1 平均軸系下穩(wěn)定性分析

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

圖4所示為轉速對飛行器穩(wěn)定性影響。由圖4(a)知，α=0.2下，考慮旋轉時，由于陀螺效應影響，正、反進動頻率(由于反進動頻率為負，這里取反進動絕對值進行分析)呈分離趨勢，1階反進動頻率與剛體模態(tài)頻率在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)區(qū)間，虛部為正數(shù)，發(fā)生動態(tài)失穩(wěn)，臨界推力降低，失穩(wěn)區(qū)域變大，引起剛性運動和彈性振動耦合振動。

當α=1.0時，由圖4(c)和圖4(d)知，1階反進動頻率與剛體模態(tài)頻率發(fā)生耦合，2階正進動頻率逐漸減小，在(0.43Pcr0，0.76Pcr0)區(qū)間，模態(tài)頻率變?yōu)榉橇悖o態(tài)失穩(wěn)變成動態(tài)失穩(wěn)，靜態(tài)失穩(wěn)推力變?yōu)閯討B(tài)失穩(wěn)推力，失穩(wěn)區(qū)域變大。

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

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

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

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

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

2.2 平均軸系和瞬態(tài)坐標系模型對比

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

表2 臨界推力

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

將平均軸系模型和瞬態(tài)坐標系模型進行仿真，結果如圖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時飛行器的穩(wěn)定性，可見瞬態(tài)坐標系模型產生剛體模態(tài)，發(fā)生靜態(tài)失穩(wěn)，臨界推力為0.31Pcr0，與平均軸系模型相比，穩(wěn)定性降低。

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

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

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

3 結論

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

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

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