潘成龍，榮吉利，徐天富，項(xiàng)大林

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

0 引 言

大長(zhǎng)徑比和高推重比自旋飛行器彈性效應(yīng)明顯，耦合問題突出，推力穩(wěn)定性問題備受關(guān)注。文獻(xiàn)[1]采用梁結(jié)構(gòu)振動(dòng)理論建立柔性彈箭運(yùn)動(dòng)模型，分析了推力作用對(duì)系統(tǒng)振動(dòng)頻率和系統(tǒng)穩(wěn)定性的影響。文獻(xiàn)[2-4]采用Timoshenko梁為模型，考慮陀螺效應(yīng)和剪切效應(yīng)，分析了隨動(dòng)推力作用下柔性自旋飛行器穩(wěn)定性問題。文獻(xiàn)[5]研究了變推力作用下系統(tǒng)質(zhì)量偏心因素對(duì)自旋飛行器系統(tǒng)動(dòng)力特性的影響。文獻(xiàn)[6]采用變質(zhì)量系統(tǒng)的方法，分析了旋轉(zhuǎn)固體火箭發(fā)動(dòng)機(jī)工作過程中的章動(dòng)不穩(wěn)定性問題。除此之外，飛行器在飛行時(shí)是一個(gè)典型的時(shí)變系統(tǒng)，燃料消耗引起變質(zhì)量效應(yīng)，對(duì)飛行過程有很大影響，其振動(dòng)規(guī)律不同于恒定質(zhì)量系統(tǒng)[7]。文獻(xiàn)[8]以Kane運(yùn)動(dòng)方程為基礎(chǔ)，建立考慮時(shí)變質(zhì)量和幾何剛度的柔性火箭結(jié)構(gòu)動(dòng)力學(xué)方程。文獻(xiàn)[9]運(yùn)用拉格朗日方程建立時(shí)變質(zhì)量柔性火箭動(dòng)力學(xué)方程，同時(shí)考慮了推進(jìn)劑質(zhì)量減少而導(dǎo)致質(zhì)心移動(dòng)的復(fù)雜情況。文獻(xiàn)[10-11]采用改進(jìn)歐拉中點(diǎn)辛差分格式、模態(tài)疊加方法和自適應(yīng)Newmark法，研究變質(zhì)量系統(tǒng)動(dòng)態(tài)的動(dòng)態(tài)響應(yīng)，證明了自適應(yīng)Newmark法對(duì)時(shí)變系統(tǒng)的適應(yīng)性。

本文在國家自然科學(xué)基金項(xiàng)目“自旋飛行器橫向振動(dòng)的動(dòng)力學(xué)與控制機(jī)理研究”工作基礎(chǔ)上，以Timoshenko梁為模型，考慮陀螺效應(yīng)和剪切效應(yīng)，基于有限元法和穩(wěn)定模態(tài)基底法[12]建立變質(zhì)量系統(tǒng)振動(dòng)方程，分析質(zhì)量變化、推力作用和自旋轉(zhuǎn)速對(duì)飛行器穩(wěn)定性的影響，采用自適應(yīng)Newmark法求解系統(tǒng)的振動(dòng)響應(yīng)。

1 動(dòng)力學(xué)方程

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

圖1 彈體模型圖

如圖1所示，建立準(zhǔn)彈體坐標(biāo)系Oxyz和隨體坐標(biāo)系O′ξηζ，忽略軸向變形，系統(tǒng)的動(dòng)能為[3]：

(1)

式中：uy，uz分別為梁截面y向，z向的橫向位移，θy，θz分別為該截面的轉(zhuǎn)角，Ω為自旋轉(zhuǎn)速，A為截面面積，ρ為軸段微元質(zhì)量密度，I為軸段微元截面的慣性矩，lb為梁長(zhǎng)。

考慮剪切影響，系統(tǒng)的彈性勢(shì)能為[3]：

(2)

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

隨動(dòng)推力做的功[3]：

(3)

軸向力PN：

(4)

在尾部，隨動(dòng)推力非保守部分做的虛功[3]：

δWP=Pθz(0,t)δuy(0,t)-Pθy(0,t)δuz(0,t)

(5)

1.2 振動(dòng)方程

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

(6)

由于燃料消耗，質(zhì)量隨時(shí)間變化，Φ和Ψ與時(shí)間有關(guān)。采用文獻(xiàn)[12]中提出的穩(wěn)定模態(tài)基底法：

(7)

將式(7)代入式(6)中：

(8)

式中：η1(t)=B(t)q1(t)，η2(t)=B(t)q2(t)。

非保守系統(tǒng)Lagrange拉格朗日方程的一般形式如下：

(9)

式中：η和Q分別為廣義坐標(biāo)和包括非保守力在內(nèi)的廣義力。取η分別為η1和η2,由此得到單元的振動(dòng)方程：

(10)

K=K1-KPN1-KP1

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

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

以長(zhǎng)徑比為25.3的等截面細(xì)長(zhǎng)回轉(zhuǎn)梁模型作為仿真模型，整體分為5段，推力作用在軸段1的左端，初始時(shí)刻相關(guān)參數(shù)見表1。圖2和3分別給出主動(dòng)段初至末時(shí)刻單位長(zhǎng)度質(zhì)量及0～3 s燃料消耗。

表1回轉(zhuǎn)梁模型參數(shù)

圖2 單位長(zhǎng)度質(zhì)量

圖3 主動(dòng)段飛行器質(zhì)量

將式(10)簡(jiǎn)化為：

(19)

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

(20)

2.1 頻率分析

由圖4知，隨著質(zhì)量的消耗，系統(tǒng)的固有頻率逐漸升高；與無推力時(shí)相比，推力能夠降低系統(tǒng)固有頻率。

圖4 固有頻率變化

圖5 不同推力作用下進(jìn)動(dòng)頻率的變化

圖5分析了質(zhì)量和恒定推力對(duì)進(jìn)動(dòng)頻率的影響。取恒定推力P=3.2×106N，Ω=102 rad·s-1，推力使正進(jìn)動(dòng)頻率減少，負(fù)反進(jìn)動(dòng)頻率減少，隨著質(zhì)量消耗，一階和二階正進(jìn)動(dòng)頻率和負(fù)反進(jìn)動(dòng)頻率都升高。在恒定推力和質(zhì)量變化共同作用下，一階進(jìn)動(dòng)頻率變化斜率越來越小，在t>2.5 s，接近零，進(jìn)動(dòng)頻率趨于不變。

圖6 推力變化

圖7 不同工況下進(jìn)動(dòng)頻率變化

圖7分析了質(zhì)量和變推力對(duì)進(jìn)動(dòng)頻率的影響，推力變化如圖6所示。取Ω=102 rad·s-1，由圖7(a)知，工況1推力增大時(shí)，一階正進(jìn)動(dòng)頻率和負(fù)反進(jìn)動(dòng)頻率先增大后減小，整體呈增大趨勢(shì)，說明這段時(shí)間質(zhì)量消耗對(duì)進(jìn)動(dòng)頻率影響大于推力作用，工況2在推力增大段，一階正進(jìn)動(dòng)頻率和負(fù)反進(jìn)動(dòng)頻率先增大后減小，整體呈減小趨勢(shì)，推力作用更明顯；工況1推力不變和減小時(shí)，一階正進(jìn)動(dòng)頻率和負(fù)反進(jìn)動(dòng)頻率逐漸正增大，推力減小時(shí)，進(jìn)動(dòng)頻率變化斜率要大于推力不變時(shí)。由圖7(b)知，二階正進(jìn)動(dòng)頻率和負(fù)反進(jìn)動(dòng)頻率逐漸增加，進(jìn)動(dòng)頻率變化斜率受推力變化影響。

圖8分析了質(zhì)量和推力對(duì)臨界轉(zhuǎn)速ωcr的影響，ωcr采用文獻(xiàn)[14]中方法計(jì)算得到。無推力時(shí)，隨著質(zhì)量消耗，臨界轉(zhuǎn)速逐漸增大；恒定推力作用時(shí)，臨界轉(zhuǎn)速變小，臨界轉(zhuǎn)速變化斜率越來越小，在t>2.5 s，接近零，臨界轉(zhuǎn)速率趨于不變；工況1和2推力增大時(shí)，臨界轉(zhuǎn)速先增大后減小，工況1推力不變和減小時(shí)，臨界轉(zhuǎn)速逐漸增大。

圖8 臨界轉(zhuǎn)速

圖9 質(zhì)心位置

圖9中，初時(shí)刻，質(zhì)心到尾部的距離為2.64 m，隨著質(zhì)量減少，質(zhì)心向頭部移動(dòng)，移動(dòng)速度變大，末時(shí)刻，到尾部的距離變?yōu)?.5 m。

2.2 振動(dòng)響應(yīng)分析

經(jīng)典Newmark法的迭代計(jì)算公式為：

(21)

(22)

在時(shí)變系統(tǒng)的振動(dòng)問題中，參數(shù)隨時(shí)間變化，采用自適應(yīng)Newmark法[10]進(jìn)行求解，步長(zhǎng)h和β為：

(23)

(24)

(25)

(26)

式中：h(t)為自適應(yīng)步長(zhǎng);κ為步長(zhǎng)控制因子；ω(t)為最高階振動(dòng)角頻率；ξ為系統(tǒng)模態(tài)阻尼比。

采用自適應(yīng)Newmark法分析由激振力引起的變質(zhì)量系統(tǒng)橫向振動(dòng)響應(yīng)問題，激振力作用于飛行器尾部y,z向，大小均為2000 N，激振頻率與自旋轉(zhuǎn)速相同。

圖10 不同轉(zhuǎn)速下尾部位移響應(yīng)

分別計(jì)算自旋轉(zhuǎn)速Ω為90 rad·s-1,101 rad·s-1,112 rad·s-1時(shí)橫向位移響應(yīng)情況。由圖8知，恒定推力P=3.2×106N時(shí)，ωcr在初和末時(shí)刻分別為95.8 rad·s-1和107.2 rad·s-1。圖10(a)中自旋轉(zhuǎn)速小于ωcr，受質(zhì)量消耗影響，尾部振幅減少；圖10(b)中振幅逐漸增大，在1.8 s時(shí)振幅達(dá)到最大，后逐漸減小，由圖8可知，在1.1 s時(shí)，自旋轉(zhuǎn)速等于臨界轉(zhuǎn)速，系統(tǒng)發(fā)生短暫共振,共振響應(yīng)峰值最大時(shí)刻與自旋轉(zhuǎn)速和臨界轉(zhuǎn)速重合時(shí)刻相比延遲0.7 s；圖10(c)中質(zhì)量消耗，臨界轉(zhuǎn)速增加，趨向于自旋轉(zhuǎn)速，導(dǎo)致尾部振幅逐漸增加。

3 結(jié) 論

以柔性自旋飛行器為研究對(duì)象，采用變質(zhì)量Timoshenko梁為模型，基于有限元法和穩(wěn)定模態(tài)基底法建立時(shí)變系統(tǒng)橫向振動(dòng)方程，通過算例仿真，分析變質(zhì)量、推力作用和轉(zhuǎn)速對(duì)飛行器穩(wěn)定性響，采用自適應(yīng)Newmark法，對(duì)系統(tǒng)橫向振動(dòng)響應(yīng)進(jìn)行分析，得到以下結(jié)論：

1)隨著質(zhì)量減少，系統(tǒng)固有頻率、進(jìn)動(dòng)頻率和臨界轉(zhuǎn)速逐漸增加，推力作用能夠減小系統(tǒng)剛度，不考慮質(zhì)量變化時(shí)，正進(jìn)動(dòng)和負(fù)反進(jìn)動(dòng)頻率減小，推力和質(zhì)量共同作用時(shí)，正進(jìn)動(dòng)和負(fù)反進(jìn)動(dòng)頻率增大，說明質(zhì)量消耗對(duì)進(jìn)動(dòng)頻率影響大于推力作用。

2)當(dāng)自旋轉(zhuǎn)速小于臨界轉(zhuǎn)速時(shí)，受質(zhì)量影響振幅減?。豢拷R界轉(zhuǎn)速時(shí)，振幅增加，等于臨界轉(zhuǎn)速頻率時(shí)發(fā)生短暫共振，共振響應(yīng)峰值最大時(shí)刻與轉(zhuǎn)速重合時(shí)刻相比稍有延遲。