齊 麟，楊 亞，陳意芬，殷 瑋，趙宏宇

（1.上海交通大學(xué)航空航天學(xué)院，上海 200240；2.上海機(jī)電工程研究所，上海 201109）

細(xì)長(zhǎng)體構(gòu)型在飛行器外形設(shè)計(jì)中有著廣泛的應(yīng)用，如運(yùn)載火箭、導(dǎo)彈及飛機(jī)機(jī)身。細(xì)長(zhǎng)體外形的自身性質(zhì)導(dǎo)致此類(lèi)飛行器在飛行過(guò)程中容易產(chǎn)生彈性形變，并且表現(xiàn)出復(fù)雜的振動(dòng)特性。振動(dòng)模態(tài)呈現(xiàn)出低頻、密頻的特點(diǎn)［1］。這些復(fù)雜密集的低頻振動(dòng)模態(tài)在控制系統(tǒng)的設(shè)計(jì)過(guò)程中起著重要的作用。由于模態(tài)頻率降低至控制系統(tǒng)帶寬內(nèi)，從而導(dǎo)致控制系統(tǒng)的輸出會(huì)受到彈性振動(dòng)的影響，伺服機(jī)構(gòu)輸出的控制力矩又導(dǎo)致氣動(dòng)力重新分布，進(jìn)一步影響細(xì)長(zhǎng)體變形。與此同時(shí)，隨著大量先進(jìn)新型材料在航空航天領(lǐng)域的廣泛應(yīng)用，在一定程度上增強(qiáng)了飛行器結(jié)構(gòu)參數(shù)的隨機(jī)性（如彈性模量、泊松比等）［2］。結(jié)構(gòu)隨機(jī)因素使得結(jié)構(gòu)動(dòng)力學(xué)穩(wěn)定性問(wèn)題呈現(xiàn)出新的特點(diǎn)［3］。

導(dǎo)致結(jié)構(gòu)失穩(wěn)的原因有很多［4］。其中一種為推力及氣動(dòng)力相互作用下導(dǎo)致的結(jié)構(gòu)失穩(wěn)［5］。Beal［6］開(kāi)展隨動(dòng)推力影響下的穩(wěn)定性問(wèn)題研究較早，對(duì)均勻梁模型展開(kāi)了常值及脈沖推力影響的穩(wěn)定性分析。Wu［7-8］將細(xì)長(zhǎng)體飛行器視為兩端自由的歐拉-伯努利梁，采用有限元方法研究了細(xì)長(zhǎng)飛行器在矢量控制推力作用下的失穩(wěn)機(jī)理，分析了結(jié)構(gòu)穩(wěn)定性與彈性模態(tài)的關(guān)系。Park［9］用均勻、兩端自由的Timoshenko梁代替細(xì)長(zhǎng)體，利用擴(kuò)展的Hamilton 原理，建立了描述橫向平面運(yùn)動(dòng)的動(dòng)力學(xué)方程，研究了轉(zhuǎn)動(dòng)慣量及剪切變形對(duì)穩(wěn)定性的影響。Trikha［10-11］考慮了剛性、縱向和橫向振動(dòng)模態(tài)的耦合，建立了分析細(xì)長(zhǎng)體飛行器在推力和空氣動(dòng)力作用下穩(wěn)定性的動(dòng)力學(xué)模型。采用時(shí)頻Fourier 譜-有限元法和H-P 有限元法對(duì)動(dòng)力學(xué)模型進(jìn)行離散，并對(duì)其發(fā)散和顫振現(xiàn)象進(jìn)行了分析。宋?。?2］推導(dǎo)了推力及空氣動(dòng)力作用下細(xì)長(zhǎng)體飛行器橫向運(yùn)動(dòng)方程及邊界條件，分析了橫向振動(dòng)模態(tài)出現(xiàn)低頻特性的原因。許赟等［13］將細(xì)長(zhǎng)體構(gòu)型簡(jiǎn)化為變質(zhì)量非均勻Euler-Bernoulli梁模型，將推力定義為隨動(dòng)力，依據(jù)Hamilton原理，建立了描述結(jié)構(gòu)變形的動(dòng)力學(xué)方程。通過(guò)有限元方法對(duì)方程進(jìn)行解算，研究了不同推力作用下結(jié)構(gòu)失穩(wěn)的特性。榮吉利等［14］將旋轉(zhuǎn)細(xì)長(zhǎng)體飛行器簡(jiǎn)化為帶有附加質(zhì)量的Timoshenko 梁模型，同時(shí)考慮旋轉(zhuǎn)導(dǎo)致的陀螺力矩及隨動(dòng)推力的影響，建立了結(jié)構(gòu)動(dòng)力學(xué)方程。通過(guò)有限元方法進(jìn)行了仿真分析，研究了旋轉(zhuǎn)速度、剪切變形、隨動(dòng)推力及附加質(zhì)量等因素對(duì)結(jié)構(gòu)穩(wěn)定性的影響。鄢雄偉等［15］建立了旋轉(zhuǎn)導(dǎo)彈在變推力影響下的耦合動(dòng)力學(xué)模型，研究了推力快速變化對(duì)彈體動(dòng)力穩(wěn)定性的影響，并給出了推力變化率設(shè)計(jì)邊界。

上述在細(xì)長(zhǎng)體構(gòu)型飛行器結(jié)構(gòu)穩(wěn)定性方面的研究，均建立在結(jié)構(gòu)參數(shù)為確定值的基礎(chǔ)上。然而，隨機(jī)性是材料的固有屬性，尤其隨著先進(jìn)材料的應(yīng)用及飛行任務(wù)需求的發(fā)展，隨機(jī)性的影響日益增強(qiáng)。因此，有必要對(duì)結(jié)構(gòu)隨機(jī)因素影響下的細(xì)長(zhǎng)體穩(wěn)定性問(wèn)題給予更多的關(guān)注。

本文將對(duì)結(jié)構(gòu)隨機(jī)因素及推力作用下的結(jié)構(gòu)穩(wěn)定性進(jìn)行研究。首先基于矩陣攝動(dòng)理論，對(duì)結(jié)構(gòu)隨機(jī)因素導(dǎo)致相應(yīng)參數(shù)矩陣發(fā)生改變時(shí)，細(xì)長(zhǎng)體動(dòng)力學(xué)模型進(jìn)行推導(dǎo)。然后，求解反映特征值變化的顯式表達(dá)式，并得到其期望與方差，形成一套隨機(jī)因素在推力作用下細(xì)長(zhǎng)體飛行器動(dòng)力學(xué)穩(wěn)定性問(wèn)題中影響的完整研究方法。同時(shí)，為驗(yàn)證方法有效性，對(duì)本文所得到的結(jié)果進(jìn)行了數(shù)值仿真驗(yàn)證。

1 推力作用下的細(xì)長(zhǎng)體飛行器動(dòng)力學(xué)建模

隨著對(duì)細(xì)長(zhǎng)體飛行器有效載荷及射程需求的增加，設(shè)計(jì)上往往要求飛行器具有較高的推重比。此時(shí)，推力對(duì)于細(xì)長(zhǎng)體彈箭動(dòng)力學(xué)穩(wěn)定性的影響需要給予考慮。首先，對(duì)推力作用下的動(dòng)力學(xué)穩(wěn)定性分析模型進(jìn)行推導(dǎo)。飛行器受力如圖1所示，其中fp為發(fā)動(dòng)機(jī)推力，fadw為氣動(dòng)阻力。

圖1 飛行器受力Fig.1 Definition of force of slender vehicle

將彈體橫向位移w(x，t)按振動(dòng)模態(tài)進(jìn)行展開(kāi)，有

式（1）中，?i(x)為振動(dòng)模態(tài)，qi(t)為廣義坐標(biāo)。不考慮結(jié)構(gòu)阻尼，有

式（2）中，及分別為質(zhì)量矩陣、結(jié)構(gòu)剛度矩陣及推力引起的附加剛度矩陣，具體形式如下：

其中，m(x)為質(zhì)量密度；E為彈性模量；I為慣性矩；N(x)為軸向力，具體形式如下：

對(duì)于動(dòng)力學(xué)穩(wěn)定性問(wèn)題，可將廣義坐標(biāo)qi(t)表示為

式（4）中，為特征值。將式（4）帶入式（3）可得

式（5）中，

由式（5）可知，A，B矩陣均為非對(duì)稱(chēng)矩陣，所以特征值為復(fù)數(shù)，形式如實(shí)部表示振動(dòng)頻率，虛部表示振動(dòng)阻尼。根據(jù)穩(wěn)定性概念，當(dāng)實(shí)部＜0 時(shí)，結(jié)構(gòu)振動(dòng)會(huì)收斂，系統(tǒng)為穩(wěn)定狀態(tài)。當(dāng)實(shí)部＞0時(shí)，結(jié)構(gòu)振動(dòng)會(huì)發(fā)散，則系統(tǒng)不穩(wěn)定。當(dāng)實(shí)部=0 時(shí)，結(jié)構(gòu)振動(dòng)位移為等幅振蕩，系統(tǒng)處于臨界狀態(tài)。此時(shí)的推力值即為臨界推力。臨界推力的計(jì)算對(duì)于結(jié)構(gòu)穩(wěn)定來(lái)說(shuō)具有重要的意義，它決定了飛行包線的形狀以及任務(wù)剖面的設(shè)計(jì)。

結(jié)構(gòu)隨機(jī)因素的隨機(jī)性會(huì)使得振動(dòng)模態(tài)?i(x)具有一定的隨機(jī)性。而振動(dòng)模態(tài)的變化會(huì)讓及發(fā)生相應(yīng)的改變，從而使得A，B矩陣具有相應(yīng)的隨機(jī)性，由此計(jì)算得到的推力臨界值也為隨機(jī)參量。

結(jié)構(gòu)隨機(jī)因素的改變會(huì)引起動(dòng)力學(xué)穩(wěn)定邊界的變化，若采用傳統(tǒng)的分析方式研究穩(wěn)定邊界的變化規(guī)律，則結(jié)構(gòu)隨機(jī)因素每變化一次就要對(duì)細(xì)長(zhǎng)體飛行器的廣義特征值問(wèn)題求解一次，需進(jìn)行多次迭代計(jì)算，分析過(guò)程非常耗時(shí)和麻煩。因此，本文采用矩陣攝動(dòng)理論建立一種能夠快速分析細(xì)長(zhǎng)體飛行器動(dòng)力學(xué)穩(wěn)定邊界變化規(guī)律的分析方法，以提高計(jì)算效率。

下面將基于矩陣攝動(dòng)理論，對(duì)考慮結(jié)構(gòu)隨機(jī)因素影響時(shí)的細(xì)長(zhǎng)體飛行器動(dòng)力學(xué)穩(wěn)定性問(wèn)題進(jìn)行建模，并求解特征值的解及其期望與方差。

2 結(jié)構(gòu)隨機(jī)因素影響下的穩(wěn)定性問(wèn)題建模

2.1 二階遞推模型建立

對(duì)于描述細(xì)長(zhǎng)體飛行器動(dòng)力學(xué)穩(wěn)定性的特征值問(wèn)題，如式（5），與其相對(duì)應(yīng)的伴隨特征值問(wèn)題可表示為

易證特征值問(wèn)題式（5）及伴隨特征值問(wèn)題式（6）具有相同的特征值。并且模態(tài)向量和滿足以下條件：

由于彈性模量、泊松比、幾何尺寸及質(zhì)量密度的結(jié)構(gòu)參數(shù)具有隨機(jī)性，因此矩陣A，B為隨機(jī)矩陣。若采用上一節(jié)的方式對(duì)A，B矩陣進(jìn)行泰勒展開(kāi)，則由于其A，B矩陣的特點(diǎn)，很難得到特征值關(guān)于隨機(jī)因素的解析解。因此，本小節(jié)采用矩陣攝動(dòng)方法研究隨機(jī)因素對(duì)特征值的影響，即研究隨機(jī)因素對(duì)細(xì)長(zhǎng)體飛行器在推力作用下動(dòng)力學(xué)穩(wěn)定性的影響。

當(dāng)結(jié)構(gòu)隨機(jī)因素產(chǎn)生有限變化時(shí)，質(zhì)量矩陣、結(jié)構(gòu)剛度矩陣及附加矩陣也會(huì)有所改變，其變化可表示為

相應(yīng)的有

式（8）-（9）中，ε表示攝動(dòng)系數(shù)，為小參數(shù)；及分別為質(zhì)量矩陣、結(jié)構(gòu)剛度矩陣及附加矩陣在隨機(jī)因素均值處的取值。分別表示隨機(jī)因素引起的變化。當(dāng)ε→0 時(shí)，即為確定性的復(fù)特征值問(wèn)題。

由于上述矩陣發(fā)生變化，因此，特征值及模態(tài)向量也會(huì)發(fā)生較小變化。將特征值及模態(tài)向量按小參數(shù)ε展開(kāi)成冪級(jí)數(shù)，有

模態(tài)向量滿足如下正則化條件：

將式（8）-（10）帶入式（5）、（6）及式（11），可建立求解一階及二階攝動(dòng)的遞推模型。將遞推公式按小參數(shù)ε的階次進(jìn)行整理。

零階方程為

一階方程為

二階方程為

通過(guò)求解一階攝動(dòng)方程可以得到特征值及模態(tài)向量的一階攝動(dòng)值。利用二階攝動(dòng)方程可以求出特征值及特征向量的二階攝動(dòng)。

2.2 一階攝動(dòng)值求解

將模態(tài)向量的一階攝動(dòng)按展開(kāi)，即

式（15）中，G0矩陣為待定系數(shù)矩陣。將式（15）帶入式（13）的第一式有

將左乘式（16），得

由模態(tài)向量的正交性可得

式（18）中

同時(shí)，待定系數(shù)矩陣G0的非對(duì)角元素可表示為

與式（15）相類(lèi)似，將模態(tài)向量按展開(kāi)，即

式（21）中，H0為待定系數(shù)矩陣。將式（21）帶入到式（13）的第二式可得

左乘上式后得

根據(jù)模態(tài)向量的正交化性質(zhì)，上式可整理為

式（24）中

對(duì)于待定系數(shù)矩陣的對(duì)角元素，可由模態(tài)向量的正交條件求得。將式（15）及式（21）帶入到式（13）的第四式，可得

整理后可得

由于=1，因此上式可改寫(xiě)為

設(shè)Gi0i=Hi0i，因此有

通過(guò)求解式（19）、（20）、（25）、（29）后可以得到特征值及模態(tài)向量的一階攝動(dòng)值。當(dāng)結(jié)構(gòu)隨機(jī)因素的變異系數(shù)較小時(shí)，即結(jié)構(gòu)參數(shù)發(fā)生較小變化，一階攝動(dòng)可以達(dá)到較理想的精度。但當(dāng)結(jié)構(gòu)隨機(jī)因素具有較大的變異系數(shù)時(shí)，需要對(duì)二階攝動(dòng)進(jìn)行求解。

2.3 二階攝動(dòng)值求解

對(duì)模態(tài)向量的二階攝動(dòng)按展開(kāi)，即

將式（30）帶入到式（14）中的第一式可得

左乘后得

根據(jù)模態(tài)向量的正交化條件，式（32）可整理為

式（33）中

從而待定系數(shù)矩陣G1的非對(duì)角元素可表示為

與模態(tài)向量的二階攝動(dòng)推導(dǎo)過(guò)程相類(lèi)似，下面推導(dǎo)模態(tài)向量的二階攝動(dòng)值。將按展開(kāi)，即

將式（36）帶入式（14）的第二式，可得

左乘式（37）后可得

利用模態(tài)向量的正交性，式（38）可整理為

式（39）中

由于為對(duì)角陣，所以待定系數(shù)矩陣H1的非對(duì)角元素可表示為

對(duì)于矩陣G1及H1的對(duì)角元素可通過(guò)模態(tài)向量的正交特性求得，將式（30）及式（36）帶入式（14）中的第四式可得

整理后可得

根據(jù)模態(tài)向量的正交性，并設(shè)Gi1i=Hi1i，可求得矩陣G1及H1的對(duì)角元素為

通過(guò)求解式（34）、（35）、（40）、（43）可得到特征值及模態(tài)向量的二階攝動(dòng)。將特征值及模態(tài)向量的一階及二階攝動(dòng)值帶入式（10）即可得到質(zhì)量及剛度矩陣發(fā)生變化時(shí)的特征值及模態(tài)向量。

前述推導(dǎo)得到了當(dāng)矩陣因結(jié)構(gòu)隨機(jī)因素而產(chǎn)生的一定的攝動(dòng)時(shí)，特征值及模態(tài)向量的變化結(jié)果。接下來(lái)將基于特征值的二階攝動(dòng)展開(kāi)，對(duì)特征值的數(shù)學(xué)期望及方差的顯式表達(dá)式進(jìn)行推導(dǎo)。

2.4 特征值的期望與方差

由于結(jié)構(gòu)隨機(jī)因素的影響，攝動(dòng)量及均可視為隨機(jī)矩陣，因此特征值也具有隨機(jī)性。為研究其統(tǒng)計(jì)特性，求取其數(shù)學(xué)期望與方差，本部分假設(shè)矩陣A0、B0、A1、B1中各元素統(tǒng)計(jì)特性已知，基于前述得到特征值的二階攝動(dòng)展開(kāi)式，對(duì)其期望及方差進(jìn)行推導(dǎo)。

依據(jù)式（10），第i個(gè)特征值的數(shù)學(xué)期望為

由于一階攝動(dòng)量及二階攝動(dòng)量的期望為0，因此，特征值的期望為

考慮特征值的二階攝動(dòng)結(jié)果，假設(shè)特征值λˉ、一階攝動(dòng)量及二階攝動(dòng)量相互獨(dú)立，則方差可表示為

假設(shè)A0、B0相互獨(dú)立，且E(A0)=0、E(B0)=0，則式（47）中展開(kāi)后的交叉項(xiàng)的期望為0，因此式（48）可寫(xiě)為

將式（49）的平方項(xiàng)展開(kāi)為標(biāo)量形式，有

將式（50）-（51）帶入式（49）后有

假設(shè)矩陣A0、B0中各元素相互獨(dú)立，則式（52）可整理為

綜合式（46）、（53），即可得到特征值的方差。

至此，二階精度的特征值期望及具有一階精度的特征值方差的計(jì)算模型推導(dǎo)完成。為對(duì)本小節(jié)推導(dǎo)的當(dāng)質(zhì)量矩陣、剛度矩陣變化時(shí)的特征值的解以及方差有效性進(jìn)行驗(yàn)證，下面將對(duì)上述推導(dǎo)進(jìn)行數(shù)值仿真分析。

3 數(shù)值仿真驗(yàn)證及分析

在本文基礎(chǔ)上可以進(jìn)一步分析推力臨界值在結(jié)構(gòu)隨機(jī)因素作用下的變化。本小節(jié)將計(jì)算當(dāng)結(jié)構(gòu)隨機(jī)因素導(dǎo)致質(zhì)量矩陣及剛度矩陣發(fā)生變化時(shí)，特征值λ的解，并與矩陣攝動(dòng)變化后的理論計(jì)算結(jié)果進(jìn)行對(duì)比，仿真條件如表1所示。

表1 仿真條件Tab.1 Simulation data of vehicle

為驗(yàn)證本文方法的有效性，當(dāng)攝動(dòng)系數(shù)分別為0.01、0.02、0.03、0.04及0.05時(shí)，進(jìn)行數(shù)值仿真計(jì)算，計(jì)算結(jié)果如表2及表3所示。表2及表3分別為一階攝動(dòng)方程及二階攝動(dòng)方程的計(jì)算結(jié)果與理論值的對(duì)比。從表中可以看出，本文所建立的方法能夠較好地計(jì)算質(zhì)量矩陣及剛度矩陣改變時(shí)，特征值的變化。并且二階攝動(dòng)得到的結(jié)果整體優(yōu)于一階攝動(dòng)。在攝動(dòng)系數(shù)較小時(shí)，即ε≤0.03時(shí)，一階攝動(dòng)所得到的結(jié)果也具有一定的精度。表4分別給出了攝動(dòng)值為0.01 時(shí)的理論值、一階攝動(dòng)值及二階攝動(dòng)值的計(jì)算時(shí)間，計(jì)算條件為CPU：Intel酷睿i5 7200U、2.5 GHz主頻，內(nèi)存：8 G。一階攝動(dòng)及二階攝動(dòng)計(jì)算時(shí)間均小于理論值，并且一階計(jì)算時(shí)間約為二階計(jì)算的1/2，約為理論值計(jì)算的1/4。由于一階攝動(dòng)方程的計(jì)算量明顯小于二階攝動(dòng)。因此，在質(zhì)量矩陣及剛度矩陣變化較小時(shí)，可采用一階攝動(dòng)方程進(jìn)行相關(guān)計(jì)算，從而節(jié)省計(jì)算時(shí)間。但當(dāng)ε＞0.03 時(shí)，一階攝動(dòng)的計(jì)算結(jié)果誤差較大，而二階攝動(dòng)依然保持較好的誤差精度。

表2 一階攝動(dòng)計(jì)算結(jié)果對(duì)比Tab.2 The calculation result comparison of the first order perturbation

表3 二階攝動(dòng)計(jì)算結(jié)果對(duì)比Tab.3 The calculation result comparison of the second order perturbation

表4 計(jì)算時(shí)間對(duì)比Tab.4 The comparison of computing time

為研究二階攝動(dòng)方程解的精度及其所能適應(yīng)的攝動(dòng)系數(shù)變化范圍，下面將進(jìn)一步改變攝動(dòng)系數(shù)，研究本方法的適用性。

當(dāng)攝動(dòng)系數(shù)ε從0 變化到0.1 時(shí)，本文計(jì)算所得的二階特征值與理論值的比較如圖2所示。圖2（a）及（b）分別為特征值實(shí)部及虛部隨變異系數(shù)的變化曲線。圖2（c）描述的是變異系數(shù)變化引起的特征值變化，橫坐標(biāo)為特征值實(shí)部，縱坐標(biāo)為虛部。從圖2可以看出，隨著攝動(dòng)系數(shù)的增大，本文計(jì)算值與理論值的誤差逐漸增大。當(dāng)ε＜0.05 時(shí)，即圖2中第一個(gè)方框中的部分，表明本文計(jì)算結(jié)果具有較好的準(zhǔn)確性。然而當(dāng)ε＞0.05 時(shí)，隨著ε的增大，本文計(jì)算結(jié)果逐漸偏離理論值。這說(shuō)明本文計(jì)算結(jié)果的準(zhǔn)確性建立在攝動(dòng)系數(shù)較小的前提上，當(dāng)ε增大到一定值時(shí)，本文方法便不再適用。通過(guò)圖2（a）、（b）可以發(fā)現(xiàn)，當(dāng)不考慮結(jié)構(gòu)阻尼時(shí)，質(zhì)量矩陣與剛度矩陣的變化對(duì)特征值虛部的影響較小。

圖2 不同ε時(shí)，本文計(jì)算結(jié)果與理論值的對(duì)比Fig.2 Comparison between the calculated results and theoretical values for different ε

為了驗(yàn)證本文推力臨界值結(jié)果的有效性，與參考文獻(xiàn)［7］中的方法所得到的結(jié)果進(jìn)行了對(duì)比，計(jì)算結(jié)果如表5及圖3。結(jié)果表明本文所得的推力臨界值具有較好的準(zhǔn)確性。至此已對(duì)結(jié)構(gòu)參數(shù)發(fā)生變化時(shí)，本文特征值計(jì)算結(jié)果及推力臨界值的正確性進(jìn)行了仿真驗(yàn)證。然而，推力影響下的細(xì)長(zhǎng)體飛行器動(dòng)力學(xué)穩(wěn)定性是最為關(guān)注的問(wèn)題。下面將通過(guò)本文方法研究質(zhì)量矩陣、剛度矩陣的變化對(duì)推力臨界值的影響。

表5 推力臨界值對(duì)比Tab.5 Comparison of critical thrust values

圖3 特征值隨推力的變化Fig.3 Variation of eigenvalue with thrust

根據(jù)前述仿真分析，當(dāng)攝動(dòng)系數(shù)ε的變化范圍為0 →0.05 時(shí)，本文特征值計(jì)算結(jié)果較為準(zhǔn)確。因此，選擇ε的攝動(dòng)范圍為0 →0.05，研究質(zhì)量矩陣、剛度矩陣發(fā)生小幅度改變時(shí)，推力臨界值的變化。

圖4表示的是攝動(dòng)系數(shù)ε的變化對(duì)推力臨界值的影響。圖4（a）為ε從0 遞增到0.05 時(shí)，特征值實(shí)部隨推力的變化?？梢钥闯霎?dāng)質(zhì)量與剛度矩陣發(fā)生變化時(shí)，推力臨界值也會(huì)產(chǎn)生相應(yīng)的改變。攝動(dòng)系數(shù)ε從0遞增到0.05，推力臨界值從23.75 kN 變化到26.13 kN，變化幅度為10%，并且兩者變化呈線性關(guān)系，如圖5所示。圖4（b）為特征值虛部的變化。可以看出，在推力臨界值附近，虛部的變化趨勢(shì)發(fā)生轉(zhuǎn)折。對(duì)比圖4（a）和（b）發(fā)現(xiàn)，質(zhì)量矩陣及剛度矩陣的變化對(duì)特征值虛部的影響小于特征值實(shí)部。這是因?yàn)樵谟?jì)算中忽略了結(jié)構(gòu)阻尼矩陣，而特征值的虛部代表了振動(dòng)特性的阻尼，因此虛部的變化幅度小于實(shí)部。

圖4 ε變化對(duì)推力臨界值的影響Fig.4 Influence of ε on the thrust critical value

圖5 ε與推力臨界值的關(guān)系Fig.5 Relationship between ε and thrust critical value

彈性模量及質(zhì)量密度等結(jié)構(gòu)參數(shù)所具有的隨機(jī)特性會(huì)帶來(lái)細(xì)長(zhǎng)體飛行器固有頻率及位移模態(tài)的變化，這些變化會(huì)導(dǎo)致動(dòng)力學(xué)穩(wěn)定性分析中的質(zhì)量矩陣及剛度矩陣產(chǎn)生相應(yīng)的改變，進(jìn)而影響動(dòng)力學(xué)發(fā)散的推力臨界值。綜合以上分析發(fā)現(xiàn)，即使當(dāng)質(zhì)量矩陣及發(fā)生小幅度改變時(shí)，推力臨界值也會(huì)產(chǎn)生較大幅度的變化（約為前者改變量的兩倍）。

當(dāng)質(zhì)量矩陣及發(fā)生隨機(jī)小幅度改變時(shí)，特征值及模態(tài)向量的值也就具有了隨機(jī)性。本文在假設(shè)矩陣A0、B0中各元素統(tǒng)計(jì)特性已知的情況下，基于矩陣攝動(dòng)理論推導(dǎo)了特征值的期望及方差。下面采用Monte-Carlo 方法對(duì)本文推導(dǎo)的特征值期望及方差進(jìn)行驗(yàn)證，樣本數(shù)量為2000，驗(yàn)證結(jié)果見(jiàn)表6。

從表6可以看出，本文方法計(jì)算得到的數(shù)學(xué)期望及方差具有一定的準(zhǔn)確性。然而方差的誤差整體上比期望的大，這是因?yàn)樵谇髮?dǎo)過(guò)程中，認(rèn)為矩陣A0、B0中的各元素相互獨(dú)立，忽略了其相關(guān)性。

表6 特征值期望及方差驗(yàn)證Tab.6 Expectation and variance verification of Eigenvalue

表6 特征值期望及方差驗(yàn)證Tab.6 Expectation and variance verification of Eigenvalue

方差/rad2·s-2推力期望/rad·s-1本文結(jié)果MCS 本文結(jié)果MCS-0.1201+1.2338i-0.2393+1.2802i-0.3055+1.3218i-0.3526+1.3596i 5 kN 10 kN 15 kN 20 kN-0.1250+1.2289i-0.2412+1.2574i-0.3066+1.3171i-0.3532+1.3549i 1.4785×10-5 9.2746×10-5 8.2458×10-5 8.0021×10-5 1.6109×10-4 9.3801×10-5 8.4660×10-5 8.1478×10-5

4 總 結(jié)

本文探究了結(jié)構(gòu)隨機(jī)因素對(duì)推力及氣動(dòng)力作用下的細(xì)長(zhǎng)體飛行器動(dòng)力學(xué)穩(wěn)定性問(wèn)題的影響?；诰仃嚁z動(dòng)理論對(duì)結(jié)構(gòu)隨機(jī)因素導(dǎo)致矩陣及發(fā)生相應(yīng)的改變時(shí)，細(xì)長(zhǎng)體飛行器動(dòng)力學(xué)模型進(jìn)行了推導(dǎo)，并給出了廣義特征值的解及其期望與方差的顯式表達(dá)式。形成了一套結(jié)構(gòu)隨機(jī)因素在推力作用下對(duì)細(xì)長(zhǎng)體飛行器動(dòng)力學(xué)穩(wěn)定性問(wèn)題中影響的完整研究方法，并對(duì)結(jié)果有效性進(jìn)行了仿真驗(yàn)證。同時(shí)，當(dāng)矩陣及發(fā)生相應(yīng)的改變時(shí)，對(duì)導(dǎo)致動(dòng)力學(xué)失穩(wěn)的推力臨界值的變化進(jìn)行了仿真分析。結(jié)果表明，廣義特征值的期望值及方差計(jì)算具有較好的精度，并且即使當(dāng)矩陣及發(fā)生小幅度改變時(shí)，推力臨界值也會(huì)產(chǎn)生較大幅度的變化。

相比于傳統(tǒng)的蒙特卡洛方法，本文所提出的方法具有不需要大量計(jì)算、計(jì)算成本低的優(yōu)勢(shì)，同時(shí)能夠?qū)逃蓄l率及模態(tài)向量的數(shù)字特征進(jìn)行直接計(jì)算。在結(jié)構(gòu)隨機(jī)因素導(dǎo)致質(zhì)量矩陣、剛度矩陣發(fā)生變化的情況下，不需要反復(fù)進(jìn)行特征值解算即可對(duì)臨界推力進(jìn)行分析，提高了計(jì)算效率。