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

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

0 引言

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

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

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

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

1 柔性彈箭動力學(xué)方程

1.1 動能、勢能和外力功

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

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

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

在準(zhǔn)彈體坐標(biāo)系下，R′和ΔC的表達(dá)式為

R′=x0+u，

(1)

(2)

(3)

則

(4)

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

在準(zhǔn)彈體坐標(biāo)系下，x0和u分別為

(5)

(6)

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

準(zhǔn)彈體坐標(biāo)系Oxyz和微元坐標(biāo)系O′ξηζ繞3個坐標(biāo)系的轉(zhuǎn)換矩陣為

(7)

在微元坐標(biāo)系下，軸段微元角速度[19]表示為

(8)

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

(9)

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

(10)

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

(11)

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

隨動推力作的功[13]為

(12)

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

(13)

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

(14)

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

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

(15)

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

1.2 振動方程

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

(16)

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

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

(17)

(18)

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

(19)

(20)

(21)

(22)

(23)

(24)

2 穩(wěn)定性分析

將(18)式簡化為

(25)

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

(26)

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

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

(27)

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

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

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

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

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

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

圖4 不同轉(zhuǎn)速和不同阻力特征值變化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ū)域變大，引起剛性運動和彈性振動耦合振動。

當(dāng)α=1.0時，由圖4(c)和圖4(d)知，1階反進(jìn)動頻率與剛體模態(tài)頻率發(fā)生耦合，2階正進(jìn)動頻率逐漸減小，在(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可見，隨著阻力的增大，臨界推力先緩慢變小，直線下降后逐漸變小。不考慮轉(zhuǎn)速時，隨阻力增大先發(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)推力；考慮轉(zhuǎn)速時，隨轉(zhuǎn)速增大，臨界推力逐漸減小，當(dāng)α>0.6時，轉(zhuǎn)速對臨界推力影響很小，轉(zhuǎn)速使剛體模態(tài)頻率與1階反進(jìn)動頻率在α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階頻率，剪切效應(yīng)對2階頻率影響更大，剪切效應(yīng)降低系統(tǒng)的穩(wěn)定性。

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

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

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

表2 臨界推力

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

將平均軸系模型和瞬態(tài)坐標(biāo)系模型進(jìn)行仿真，結(jié)果如圖8～圖10所示。

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

圖9 阻力作用下不同坐標(biāo)系模型特征值變化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所示為轉(zhuǎn)速和阻力為0時飛行器的穩(wěn)定性，可見瞬態(tài)坐標(biāo)系模型產(chǎn)生剛體模態(tài)，發(fā)生靜態(tài)失穩(wěn)，臨界推力為0.31Pcr0，與平均軸系模型相比，穩(wěn)定性降低。

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

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

由圖8～圖10可知：阻力較小時，瞬態(tài)坐標(biāo)系模型穩(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 結(jié)論

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

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

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