李海東，沈奕哲，侯凱宇，夏 鵬，劉陸廣，管耀耀

（1.上海機(jī)電工程研究所，上海 201109；2.上海航天技術(shù)研究院，上海 201109；3.上海航天精密機(jī)械研究所，上海 201600）

0 引 言

新一代飛行器向超高速、大攻角、高機(jī)動(dòng)、大長(zhǎng)細(xì)比和輕量化方向發(fā)展，其結(jié)構(gòu)特性發(fā)生顯著改變，更加柔性，受到擾動(dòng)后更易產(chǎn)生彈性變形［1］。傳統(tǒng)飛行器各系統(tǒng)孤立設(shè)計(jì)、分段驗(yàn)證的方法已經(jīng)不能適應(yīng)其發(fā)展的需要［2］。尤其是飛行器采用的固體火箭發(fā)動(dòng)機(jī)的輸出推力矢量方向始終與結(jié)構(gòu)的尾端相切，飛行過程中隨著細(xì)長(zhǎng)體飛行器振動(dòng)變形，軸線發(fā)生變形，推力方向也隨著振動(dòng)改變。這種隨著結(jié)構(gòu)振動(dòng)改變作用方向的力稱為隨動(dòng)推力。

早在20 世紀(jì)60 年代，細(xì)長(zhǎng)體導(dǎo)彈的彈性效應(yīng)問題就開始受到學(xué)者關(guān)注，更是當(dāng)前航天工程中的熱點(diǎn)問題。Beal 等［3］在1965 年研究理想化為均勻自由梁的彈性導(dǎo)彈在隨動(dòng)推力影響下的動(dòng)力學(xué)穩(wěn)定性問題，并考慮脈動(dòng)推力和簡(jiǎn)單的方向反饋控制系統(tǒng)的影響。Mladenov 等［4］研究鉸接自由梁結(jié)構(gòu)的振動(dòng)穩(wěn)定性問題。Joshi 等［5］研究大推力對(duì)彈體結(jié)構(gòu)穩(wěn)定性的影響。宋健［6］研究空氣阻力和推力同時(shí)作用時(shí)，火箭橫向振動(dòng)的模態(tài)特性，并對(duì)火箭模態(tài)頻率降低的原因進(jìn)行解釋。許赟等［7］將彈箭飛行器簡(jiǎn)化為非均勻梁模型，采用有限元法分析隨動(dòng)力作用下彈箭結(jié)構(gòu)的模態(tài)特性和穩(wěn)定性，研究彈箭結(jié)構(gòu)的模態(tài)特性（模態(tài)頻率和模態(tài)振型）及穩(wěn)定性受隨動(dòng)力的影響。羅夢(mèng)翔等［8］假設(shè)推力阻力平衡，研究推力以及軸向運(yùn)動(dòng)共同作用下彈體結(jié)構(gòu)的穩(wěn)定性。Esbati 等［9］研究三級(jí)運(yùn)載火箭在隨動(dòng)推力作用下的模態(tài)特性以及穩(wěn)定性。Ahmadia 等［10］對(duì)隨動(dòng)推力作用下兩端自由的鉸接梁動(dòng)力學(xué)穩(wěn)定性的變化展開研究。滕兆春等［11］建立了梁在軸向載荷作用下過屈曲橫向自由振動(dòng)的精確模型，獲得線性振動(dòng)的響應(yīng)，發(fā)現(xiàn)軸向力對(duì)過屈曲前后梁的各階固有頻率均有影響。Pourtakdoust 等［12］研究推力對(duì)柔性導(dǎo)彈彎曲特性的影響，發(fā)現(xiàn)質(zhì)量減小增加導(dǎo)彈振動(dòng)頻率，推力作用減小振動(dòng)頻率。

關(guān)于隨動(dòng)推力是否可以作為隨動(dòng)載荷加載在結(jié)構(gòu)上，并影響結(jié)構(gòu)的振動(dòng)特性及動(dòng)力學(xué)穩(wěn)定性一直存在爭(zhēng)議［1，13］。Sugiyama 等［14-16］開展固體火箭發(fā)動(dòng)機(jī)推力作用于懸臂梁自由端切線的相關(guān)試驗(yàn)，驗(yàn)證了工程應(yīng)用中考慮隨動(dòng)推力的必要性。但是關(guān)于隨動(dòng)推力是怎么影響導(dǎo)彈細(xì)長(zhǎng)體結(jié)構(gòu)模態(tài)特性和穩(wěn)定性的具體問題，前文學(xué)者多數(shù)只考慮隨動(dòng)推力軸向分量的影響，未考慮隨動(dòng)推力橫向分量的影響，部分引入隨動(dòng)推力橫向分量?jī)H進(jìn)行定性分析，對(duì)橫向分量對(duì)導(dǎo)彈模態(tài)特性和穩(wěn)定性的影響未進(jìn)行定量對(duì)比分析。

本文將開展隨動(dòng)推力作用下細(xì)長(zhǎng)體結(jié)構(gòu)振動(dòng)特性數(shù)值仿真研究，對(duì)引入隨動(dòng)推力橫向分量后細(xì)長(zhǎng)體結(jié)構(gòu)模態(tài)特性、穩(wěn)定性進(jìn)行分析，并與僅考慮隨動(dòng)推力軸向分量作用時(shí)的情況進(jìn)行對(duì)比，定量分析有無(wú)隨動(dòng)推力橫向分量對(duì)模態(tài)特性的影響。

1 結(jié)構(gòu)動(dòng)力學(xué)模型

隨動(dòng)推力作用下的細(xì)長(zhǎng)體結(jié)構(gòu)示意如圖1 所示。圖中：F為隨動(dòng)推力；Ff為隨動(dòng)推力橫向分量；FN為隨動(dòng)推力軸向分量；w為橫向位移；x為細(xì)長(zhǎng)體結(jié)構(gòu)軸向以及其運(yùn)動(dòng)方向坐標(biāo)。將該細(xì)長(zhǎng)體結(jié)構(gòu)劃分為n個(gè)梁?jiǎn)卧?，考察其中長(zhǎng)為l、密度為ρ、截面積為A、抗彎剛度為EI的第e個(gè)梁?jiǎn)卧?/p>

圖1 隨動(dòng)推力作用下細(xì)長(zhǎng)體Fig.1 The diagram of the slender body under follower thrust

該梁?jiǎn)卧獎(jiǎng)幽躎表達(dá)式為

設(shè)梁?jiǎn)卧獌啥斯?jié)點(diǎn)分別為i、j，梁?jiǎn)卧?jié)點(diǎn)位移w(e)是兩端節(jié)點(diǎn)的撓度和轉(zhuǎn)角，即

式中：wi、wj為梁?jiǎn)卧獌啥斯?jié)點(diǎn)的橫向位移；θi、θj為梁?jiǎn)卧獌啥斯?jié)點(diǎn)的轉(zhuǎn)角，θi、θj為小量。

梁?jiǎn)卧臋M向位移為

將式（3）代入式（1），得

式中：M(e)為梁?jiǎn)卧|(zhì)量矩陣。

由于引入隨動(dòng)推力，該梁?jiǎn)卧獎(jiǎng)菽躒表示為

式中，P(x)為隨動(dòng)推力軸向分量在梁軸向的分布函數(shù)。

將式（3）、式（8）代入式（7），得

式中：K(e)s為梁?jiǎn)卧Y(jié)構(gòu)剛度矩陣；K(e)p為梁?jiǎn)卧鯌?yīng)力矩陣。

由式（11）可見，梁?jiǎn)卧某鯌?yīng)力矩陣由隨動(dòng)推力軸向分量引入。關(guān)于隨動(dòng)推力的影響，均只考慮了隨動(dòng)推力的軸向分量FN，但隨動(dòng)推力在橫向也存在分量Ff，如圖2所示。

圖2 隨動(dòng)推力分解Fig.2 The decomposition of follower thrust

由圖2可得

Ff作用在最后一個(gè)單元的j節(jié)點(diǎn)上，因此，廣義力f(e)為

式中，K(e)f為隨動(dòng)推力橫向分量引入的剛度矩陣。

應(yīng)用拉格朗日方程建立振動(dòng)方程為

式中：qj為廣義坐標(biāo)；f(e)j為對(duì)應(yīng)于廣義坐標(biāo)的廣義力函數(shù)。

將式（1）、式（7）帶入式（16）得

由式（13）、式（16）、式（17）可得隨動(dòng)推力作用下梁?jiǎn)卧臋M向運(yùn)動(dòng)方程為

整理得

式（20）中，隨動(dòng)推力橫向分量作為剛度項(xiàng)引入。根據(jù)各個(gè)單元之間的變形協(xié)調(diào)條件和平衡條件，得到整個(gè)細(xì)長(zhǎng)體結(jié)構(gòu)的質(zhì)量矩陣M，剛度矩陣K，從而得到整個(gè)細(xì)長(zhǎng)體結(jié)構(gòu)的橫向振動(dòng)方程為

解式（21）的廣義特征值，得到隨動(dòng)推力作用下細(xì)長(zhǎng)體結(jié)構(gòu)的模態(tài)特性。

2 模態(tài)特性及穩(wěn)定性研究

對(duì)于不同的發(fā)動(dòng)機(jī)隨動(dòng)推力F，求解其作用下細(xì)長(zhǎng)體結(jié)構(gòu)的橫向振動(dòng)方程，得到一系列廣義特征值。如果某個(gè)特征值出現(xiàn)了正實(shí)部，且此時(shí)某階特征頻率f降為0 或者某兩階頻率f1、f2重合，稱此時(shí)的隨動(dòng)推力為臨界推力，記為Fcr。

設(shè)計(jì)一個(gè)1 200 mm×30 mm×9 mm 的細(xì)長(zhǎng)體鋁桿，對(duì)隨動(dòng)推力F作用下該細(xì)長(zhǎng)體結(jié)構(gòu)穩(wěn)定性進(jìn)行分析，結(jié)果見表1。由表1 中可見，不考慮隨動(dòng)推力橫向分量情況下，細(xì)長(zhǎng)體鋁棒的不穩(wěn)定屈曲問題是經(jīng)典的壓桿穩(wěn)定問題，此時(shí)細(xì)長(zhǎng)體系統(tǒng)的屈曲臨界軸向力為180 N；考慮隨動(dòng)推力橫向分量情況下，細(xì)長(zhǎng)體鋁棒的不穩(wěn)定問題是動(dòng)力學(xué)失穩(wěn)問題，此時(shí)細(xì)長(zhǎng)體系統(tǒng)的臨界推力Fcr為1 620 N。在隨動(dòng)推力橫向分量作用下，該細(xì)長(zhǎng)體系統(tǒng)的動(dòng)力學(xué)失穩(wěn)臨界值明顯高于靜態(tài)屈曲的臨界值。

表1 細(xì)長(zhǎng)體系統(tǒng)不同情況下的臨界力值Tab.1 The critical force of the slender body structure

該細(xì)長(zhǎng)體系統(tǒng)不考慮隨動(dòng)推力橫向分量、考慮隨動(dòng)推力橫向分量情況下的頻率與推力關(guān)系如圖3和圖4 所示。不考慮隨動(dòng)推力橫向分量情況的屈曲問題中，推力作用方向始終沿著未變形前的細(xì)長(zhǎng)體軸線。由圖3 可見，這種情況下模態(tài)頻率f1、f2只會(huì)依次單調(diào)下降至0，不會(huì)出現(xiàn)系統(tǒng)頻率相互接近、模態(tài)耦合的現(xiàn)象。由圖4 可見，在考慮隨動(dòng)推力橫向分量情況下，隨著隨動(dòng)推力逐漸增大，細(xì)長(zhǎng)體系統(tǒng)1階頻率f1逐漸上升，2 階頻率f2逐漸下降，兩者呈相互靠近的趨勢(shì)，在臨界點(diǎn)上1、2 階頻率重合，證實(shí)隨著細(xì)長(zhǎng)體彎曲振動(dòng)，隨動(dòng)推力作為隨動(dòng)載荷始終作用于細(xì)長(zhǎng)體自由端切線方向，其橫向分量引起了模態(tài)間的剛度耦合。

圖3 細(xì)長(zhǎng)體系統(tǒng)頻率與推力關(guān)系曲線（不考慮隨動(dòng)推力橫向分量）Fig.3 Frequency-thrust curve of the slender body （no thrust transverse component）

圖4 細(xì)長(zhǎng)體系統(tǒng)頻率曲線（考慮隨動(dòng)推力橫向分量）Fig.4 Frequency-thrust curve of slender body （considering thrust transverse component）

3 結(jié)束語(yǔ)

1） 隨動(dòng)推力軸向分量為保守力，引起細(xì)長(zhǎng)體軸向剛度的削弱；細(xì)長(zhǎng)體彎曲振動(dòng)引起的隨動(dòng)推力橫向分量為非保守力，引起系統(tǒng)模態(tài)間的耦合，容易產(chǎn)生動(dòng)力學(xué)失穩(wěn)。

2） 對(duì)于一個(gè)給定的細(xì)長(zhǎng)體結(jié)構(gòu)，在隨動(dòng)推力作用下，系統(tǒng)的動(dòng)力學(xué)失穩(wěn)臨界值高于靜態(tài)屈曲的臨界值；在動(dòng)力學(xué)失穩(wěn)的臨界點(diǎn)上，細(xì)長(zhǎng)體系統(tǒng)的1、2階頻率重頻。