結構隨機細長體飛行器動力學建模與分析

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

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

細長體構型在飛行器外形設計中有著廣泛的應用，如運載火箭、導彈及飛機機身。細長體外形的自身性質導致此類飛行器在飛行過程中容易產生彈性形變，并且表現出復雜的振動特性。振動模態呈現出低頻、密頻的特點［1］。這些復雜密集的低頻振動模態在控制系統的設計過程中起著重要的作用。由于模態頻率降低至控制系統帶寬內，從而導致控制系統的輸出會受到彈性振動的影響，伺服機構輸出的控制力矩又導致氣動力重新分布，進一步影響細長體變形。與此同時，隨著大量先進新型材料在航空航天領域的廣泛應用，在一定程度上增強了飛行器結構參數的隨機性（如彈性模量、泊松比等）［2］。結構隨機因素使得結構動力學穩定性問題呈現出新的特點［3］。

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

上述在細長體構型飛行器結構穩定性方面的研究，均建立在結構參數為確定值的基礎上。然而，隨機性是材料的固有屬性，尤其隨著先進材料的應用及飛行任務需求的發展，隨機性的影響日益增強。因此，有必要對結構隨機因素影響下的細長體穩定性問題給予更多的關注。

本文將對結構隨機因素及推力作用下的結構穩定性進行研究。首先基于矩陣攝動理論，對結構隨機因素導致相應參數矩陣發生改變時，細長體動力學模型進行推導。然后，求解反映特征值變化的顯式表達式，并得到其期望與方差，形成一套隨機因素在推力作用下細長體飛行器動力學穩定性問題中影響的完整研究方法。同時，為驗證方法有效性，對本文所得到的結果進行了數值仿真驗證。

1 推力作用下的細長體飛行器動力學建模

隨著對細長體飛行器有效載荷及射程需求的增加，設計上往往要求飛行器具有較高的推重比。此時，推力對于細長體彈箭動力學穩定性的影響需要給予考慮。首先，對推力作用下的動力學穩定性分析模型進行推導。飛行器受力如圖1所示，其中fp為發動機推力，fadw為氣動阻力。

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

將彈體橫向位移w(x，t)按振動模態進行展開，有

式（1）中，?i(x)為振動模態，qi(t)為廣義坐標。不考慮結構阻尼，有

式（2）中，及分別為質量矩陣、結構剛度矩陣及推力引起的附加剛度矩陣，具體形式如下：

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

對于動力學穩定性問題，可將廣義坐標qi(t)表示為

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

式（5）中，

由式（5）可知，A，B矩陣均為非對稱矩陣，所以特征值為復數，形式如實部表示振動頻率，虛部表示振動阻尼。根據穩定性概念，當實部＜0 時，結構振動會收斂，系統為穩定狀態。當實部＞0時，結構振動會發散，則系統不穩定。當實部=0 時，結構振動位移為等幅振蕩，系統處于臨界狀態。此時的推力值即為臨界推力。臨界推力的計算對于結構穩定來說具有重要的意義，它決定了飛行包線的形狀以及任務剖面的設計。

結構隨機因素的隨機性會使得振動模態?i(x)具有一定的隨機性。而振動模態的變化會讓及發生相應的改變，從而使得A，B矩陣具有相應的隨機性，由此計算得到的推力臨界值也為隨機參量。

結構隨機因素的改變會引起動力學穩定邊界的變化，若采用傳統的分析方式研究穩定邊界的變化規律，則結構隨機因素每變化一次就要對細長體飛行器的廣義特征值問題求解一次，需進行多次迭代計算，分析過程非常耗時和麻煩。因此，本文采用矩陣攝動理論建立一種能夠快速分析細長體飛行器動力學穩定邊界變化規律的分析方法，以提高計算效率。

下面將基于矩陣攝動理論，對考慮結構隨機因素影響時的細長體飛行器動力學穩定性問題進行建模，并求解特征值的解及其期望與方差。

2 結構隨機因素影響下的穩定性問題建模

2.1 二階遞推模型建立

對于描述細長體飛行器動力學穩定性的特征值問題，如式（5），與其相對應的伴隨特征值問題可表示為

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

由于彈性模量、泊松比、幾何尺寸及質量密度的結構參數具有隨機性，因此矩陣A，B為隨機矩陣。若采用上一節的方式對A，B矩陣進行泰勒展開，則由于其A，B矩陣的特點，很難得到特征值關于隨機因素的解析解。因此，本小節采用矩陣攝動方法研究隨機因素對特征值的影響，即研究隨機因素對細長體飛行器在推力作用下動力學穩定性的影響。

當結構隨機因素產生有限變化時，質量矩陣、結構剛度矩陣及附加矩陣也會有所改變，其變化可表示為

相應的有

式（8）-（9）中，ε表示攝動系數，為小參數；及分別為質量矩陣、結構剛度矩陣及附加矩陣在隨機因素均值處的取值。分別表示隨機因素引起的變化。當ε→0 時，即為確定性的復特征值問題。

由于上述矩陣發生變化，因此，特征值及模態向量也會發生較小變化。將特征值及模態向量按小參數ε展開成冪級數，有

模態向量滿足如下正則化條件：

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

零階方程為

一階方程為

二階方程為

通過求解一階攝動方程可以得到特征值及模態向量的一階攝動值。利用二階攝動方程可以求出特征值及特征向量的二階攝動。

2.2 一階攝動值求解

將模態向量的一階攝動按展開，即

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

將左乘式（16），得

由模態向量的正交性可得

式（18）中

同時，待定系數矩陣G0的非對角元素可表示為

與式（15）相類似，將模態向量按展開，即

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

左乘上式后得

根據模態向量的正交化性質，上式可整理為

式（24）中

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

整理后可得

由于=1，因此上式可改寫為

設Gi0i=Hi0i，因此有

通過求解式（19）、（20）、（25）、（29）后可以得到特征值及模態向量的一階攝動值。當結構隨機因素的變異系數較小時，即結構參數發生較小變化，一階攝動可以達到較理想的精度。但當結構隨機因素具有較大的變異系數時，需要對二階攝動進行求解。

2.3 二階攝動值求解

對模態向量的二階攝動按展開，即

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

左乘后得

根據模態向量的正交化條件，式（32）可整理為

式（33）中

從而待定系數矩陣G1的非對角元素可表示為

與模態向量的二階攝動推導過程相類似，下面推導模態向量的二階攝動值。將按展開，即

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

左乘式（37）后可得

利用模態向量的正交性，式（38）可整理為

式（39）中

由于為對角陣，所以待定系數矩陣H1的非對角元素可表示為

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

整理后可得

根據模態向量的正交性，并設Gi1i=Hi1i，可求得矩陣G1及H1的對角元素為

通過求解式（34）、（35）、（40）、（43）可得到特征值及模態向量的二階攝動。將特征值及模態向量的一階及二階攝動值帶入式（10）即可得到質量及剛度矩陣發生變化時的特征值及模態向量。

前述推導得到了當矩陣因結構隨機因素而產生的一定的攝動時，特征值及模態向量的變化結果。接下來將基于特征值的二階攝動展開，對特征值的數學期望及方差的顯式表達式進行推導。

2.4 特征值的期望與方差

由于結構隨機因素的影響，攝動量及均可視為隨機矩陣，因此特征值也具有隨機性。為研究其統計特性，求取其數學期望與方差，本部分假設矩陣A0、B0、A1、B1中各元素統計特性已知，基于前述得到特征值的二階攝動展開式，對其期望及方差進行推導。

依據式（10），第i個特征值的數學期望為

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

考慮特征值的二階攝動結果，假設特征值λˉ、一階攝動量及二階攝動量相互獨立，則方差可表示為

假設A0、B0相互獨立，且E(A0)=0、E(B0)=0，則式（47）中展開后的交叉項的期望為0，因此式（48）可寫為

將式（49）的平方項展開為標量形式，有

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

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

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

至此，二階精度的特征值期望及具有一階精度的特征值方差的計算模型推導完成。為對本小節推導的當質量矩陣、剛度矩陣變化時的特征值的解以及方差有效性進行驗證，下面將對上述推導進行數值仿真分析。

3 數值仿真驗證及分析

在本文基礎上可以進一步分析推力臨界值在結構隨機因素作用下的變化。本小節將計算當結構隨機因素導致質量矩陣及剛度矩陣發生變化時，特征值λ的解，并與矩陣攝動變化后的理論計算結果進行對比，仿真條件如表1所示。

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

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

表2 一階攝動計算結果對比Tab.2 The calculation result comparison of the first order perturbation

表3 二階攝動計算結果對比Tab.3 The calculation result comparison of the second order perturbation

表4 計算時間對比Tab.4 The comparison of computing time

為研究二階攝動方程解的精度及其所能適應的攝動系數變化范圍，下面將進一步改變攝動系數，研究本方法的適用性。

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

圖2 不同ε時，本文計算結果與理論值的對比Fig.2 Comparison between the calculated results and theoretical values for different ε

為了驗證本文推力臨界值結果的有效性，與參考文獻［7］中的方法所得到的結果進行了對比，計算結果如表5及圖3。結果表明本文所得的推力臨界值具有較好的準確性。至此已對結構參數發生變化時，本文特征值計算結果及推力臨界值的正確性進行了仿真驗證。然而，推力影響下的細長體飛行器動力學穩定性是最為關注的問題。下面將通過本文方法研究質量矩陣、剛度矩陣的變化對推力臨界值的影響。

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

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

根據前述仿真分析，當攝動系數ε的變化范圍為0 →0.05 時，本文特征值計算結果較為準確。因此，選擇ε的攝動范圍為0 →0.05，研究質量矩陣、剛度矩陣發生小幅度改變時，推力臨界值的變化。

圖4表示的是攝動系數ε的變化對推力臨界值的影響。圖4（a）為ε從0 遞增到0.05 時，特征值實部隨推力的變化。可以看出當質量與剛度矩陣發生變化時，推力臨界值也會產生相應的改變。攝動系數ε從0遞增到0.05，推力臨界值從23.75 kN 變化到26.13 kN，變化幅度為10%，并且兩者變化呈線性關系，如圖5所示。圖4（b）為特征值虛部的變化。可以看出，在推力臨界值附近，虛部的變化趨勢發生轉折。對比圖4（a）和（b）發現，質量矩陣及剛度矩陣的變化對特征值虛部的影響小于特征值實部。這是因為在計算中忽略了結構阻尼矩陣，而特征值的虛部代表了振動特性的阻尼，因此虛部的變化幅度小于實部。

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

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

彈性模量及質量密度等結構參數所具有的隨機特性會帶來細長體飛行器固有頻率及位移模態的變化，這些變化會導致動力學穩定性分析中的質量矩陣及剛度矩陣產生相應的改變，進而影響動力學發散的推力臨界值。綜合以上分析發現，即使當質量矩陣及發生小幅度改變時，推力臨界值也會產生較大幅度的變化（約為前者改變量的兩倍）。

當質量矩陣及發生隨機小幅度改變時，特征值及模態向量的值也就具有了隨機性。本文在假設矩陣A0、B0中各元素統計特性已知的情況下，基于矩陣攝動理論推導了特征值的期望及方差。下面采用Monte-Carlo 方法對本文推導的特征值期望及方差進行驗證，樣本數量為2000，驗證結果見表6。

從表6可以看出，本文方法計算得到的數學期望及方差具有一定的準確性。然而方差的誤差整體上比期望的大，這是因為在求導過程中，認為矩陣A0、B0中的各元素相互獨立，忽略了其相關性。

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

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

方差/rad2·s-2推力期望/rad·s-1本文結果MCS 本文結果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 總 結

本文探究了結構隨機因素對推力及氣動力作用下的細長體飛行器動力學穩定性問題的影響。基于矩陣攝動理論對結構隨機因素導致矩陣及發生相應的改變時，細長體飛行器動力學模型進行了推導，并給出了廣義特征值的解及其期望與方差的顯式表達式。形成了一套結構隨機因素在推力作用下對細長體飛行器動力學穩定性問題中影響的完整研究方法，并對結果有效性進行了仿真驗證。同時，當矩陣及發生相應的改變時，對導致動力學失穩的推力臨界值的變化進行了仿真分析。結果表明，廣義特征值的期望值及方差計算具有較好的精度，并且即使當矩陣及發生小幅度改變時，推力臨界值也會產生較大幅度的變化。

相比于傳統的蒙特卡洛方法，本文所提出的方法具有不需要大量計算、計算成本低的優勢，同時能夠對固有頻率及模態向量的數字特征進行直接計算。在結構隨機因素導致質量矩陣、剛度矩陣發生變化的情況下，不需要反復進行特征值解算即可對臨界推力進行分析，提高了計算效率。