惠俊鵬，楊 超

（1. 中國運(yùn)載火箭技術(shù)研究院研究發(fā)展中心，北京，100076；2. 北京航空航天大學(xué)，北京，100191）

跨大氣層飛行器氣動伺服彈性穩(wěn)定性分析方法研究

惠俊鵬1,2，楊 超2

（1. 中國運(yùn)載火箭技術(shù)研究院研究發(fā)展中心，北京，100076；2. 北京航空航天大學(xué)，北京，100191）

跨大氣層飛行器通常采用翼身組合體外形，機(jī)體結(jié)構(gòu)柔性大，飛行控制系統(tǒng)通道之間交聯(lián)耦合且通頻帶寬、再入擾動因素復(fù)雜、氣動加熱嚴(yán)重，上述因素可能導(dǎo)致氣動伺服彈性或熱氣動伺服彈性問題?？紤]外界干擾不確定性、飛行器模型攝動和控制通道耦合等因素，跨大氣層飛行器需進(jìn)行多通道交聯(lián)耦合氣動伺服彈性魯棒穩(wěn)定性分析。通過彈性飛行器動力學(xué)建模、非定常氣動力擬合、伺服系統(tǒng)和飛行控制系統(tǒng)建模，建立了氣動伺服彈性系統(tǒng)的閉環(huán)模型。在此基礎(chǔ)上對 Nyquist方法、最小奇異值法以及結(jié)構(gòu)奇異值μ方法等氣動伺服彈性穩(wěn)定性分析方法進(jìn)行了分析與討論，得出相關(guān)結(jié)論。

跨大氣層飛行器；氣動伺服彈性；穩(wěn)定性；結(jié)構(gòu)奇異值；μ方法

0 引 言

飛行器氣動伺服彈性（Aeroservo-elasticity，ASE）問題是一門涉及彈性飛行器結(jié)構(gòu)、非定常氣動力以及飛行控制系統(tǒng)三者相互作用的多學(xué)科技術(shù)問題。該問題實(shí)質(zhì)為：飛行器在飛行中受擾動而發(fā)生運(yùn)動，這種運(yùn)動包含了飛行器的剛體運(yùn)動和機(jī)體結(jié)構(gòu)的彈性振動。傳感器同時(shí)接收剛體運(yùn)動信號和結(jié)構(gòu)彈性振動信號，如伺服控制系統(tǒng)通頻帶覆蓋飛行器主要固有頻率，則舵面的反饋偏轉(zhuǎn)中將額外增加一部分頻率較高的偏轉(zhuǎn)運(yùn)動（相應(yīng)的會增加一部分舵面非定常氣動力和慣性力）。這是由于控制系統(tǒng)的反饋而額外作用于飛行器的控制力，這部分控制力將影響飛行器的運(yùn)動，由此形成飛行器氣動伺服彈性系統(tǒng)，如圖1所示[1,2]。

圖1 氣動伺服彈性系統(tǒng)

氣動伺服彈性穩(wěn)定性分析是指對于給定的飛行器，研究氣動伺服彈性閉環(huán)系統(tǒng)的穩(wěn)定性。

1 跨大氣層飛行器ASE分析的特點(diǎn)

相比于其它飛行器，跨大氣層飛行器具有以下特點(diǎn)[3,4]：

a）飛行器機(jī)體尺寸大，機(jī)體結(jié)構(gòu)固有頻率較低，容易發(fā)生氣動彈性/氣動伺服彈性問題；

b）再入飛行中氣動加熱效應(yīng)顯著，可能引起飛行器結(jié)構(gòu)剛度發(fā)生變化，給飛行器氣動彈性特性帶來不利影響；

c）飛行器飛行速域、空域范圍大，飛行過程中高度及動壓變化范圍大，要求飛行器具有一定的氣動伺服彈性魯棒穩(wěn)定性；

d）飛行器通常采用翼身組合體或升力體氣動布局，氣動特性復(fù)雜，擾動因素復(fù)雜；

e）飛行器常采用多個氣動操縱面進(jìn)行三通道飛行控制，各通道之間存在相互影響作用，難以完全解耦，是典型的多通道交聯(lián)耦合系統(tǒng)；

f）飛行過程中燃料的消耗易導(dǎo)致飛行器質(zhì)量特性發(fā)生變化，給飛行器的動力學(xué)建模引入顯著攝動。

由于上述因素，跨大氣層飛行器氣動伺服彈性分析面臨如下問題：

a）飛行包線寬，飛行過程中外界擾動因素復(fù)雜；

b）受氣動加熱和燃料消耗影響，飛行器模型可能發(fā)生攝動；

c）飛行器控制通道之間存在交聯(lián)耦合，這種交聯(lián)耦合作用對系統(tǒng)穩(wěn)定性的影響往往是不利的。

因此對于跨大氣層飛行器而言，多通道交聯(lián)耦合氣動伺服彈性系統(tǒng)建模與分析技術(shù)是研制過程中需解決的重要問題。

2 跨大氣層飛行器氣動伺服彈性系統(tǒng)建模

跨大氣層飛行器氣動伺服彈性系統(tǒng)是由飛行器彈性機(jī)體子系統(tǒng)、伺服作動子系統(tǒng)、飛行控制子系統(tǒng)等組成的閉環(huán)系統(tǒng)，因此建立飛行器氣動伺服系統(tǒng)模型需要建立各個子系統(tǒng)的數(shù)學(xué)模型。

2.1 彈性機(jī)體動力學(xué)模型

綜合考慮飛行器的剛體和彈性自由度，建立彈性飛行器運(yùn)動方程。為避免采用結(jié)構(gòu)體軸系或風(fēng)軸系而出現(xiàn)剛體運(yùn)動與彈性運(yùn)動之間的慣性耦合項(xiàng)，選用一個平均軸系，可以滿足在任何時(shí)刻彈性變形引起的相對線動量和角動量均為0的條件，即[3,5,6]：

式中 d為彈性結(jié)構(gòu)上某一質(zhì)量塊dm相對于平均軸系的位置向量；δ δt為相對于平均軸系的時(shí)間導(dǎo)數(shù)。

設(shè)ijk，， 為平均軸系的單位向量，考慮彈性模態(tài)之間的正交性，可得系統(tǒng)的動能表達(dá)式為

式中 m為飛行器總質(zhì)量；Mi為第i階彈性模態(tài)的廣義質(zhì)量； ξi為廣義彈性模態(tài)坐標(biāo)；u，v，w及 ωT分別為飛行器相對于慣性系的線速度和角速度用在平均軸上的分量，其中， ωT=[ p qr ]；I為轉(zhuǎn)動慣性張量。

考慮系統(tǒng)的彈性應(yīng)變勢能，利用結(jié)構(gòu)的自由振動彈性模態(tài)，系統(tǒng)應(yīng)變能為

式中 ωi為自由振動的第i階彈性模態(tài)的固有頻率。

設(shè)Fx，F(xiàn)y，F(xiàn)z和 Mx，My，Mz分別為外力、外力矩在平均軸上的分量，應(yīng)用拉格朗日方程可得彈性飛行器的運(yùn)動方程為

式中 Qξi為對應(yīng)于第i階彈性模態(tài)的廣義力。

由于飛行器結(jié)構(gòu)彈性模態(tài)與剛體模態(tài)正交，在小擾動假設(shè)下，飛行器的彈性運(yùn)動方程為[3,6]

式中 K為廣義剛度矩陣；V0為遠(yuǎn)前方來流速度；g為重力加速度；θ0為平衡狀態(tài)處的俯仰角；M為廣義質(zhì)量矩陣，M=diag（M1，M2，…，Mn）；A為廣義非定常氣動力影響系數(shù)矩陣；φ和θ為慣性系到平均軸系的轉(zhuǎn)動歐拉角；ξ，η和δ分別為彈性模態(tài)、剛體運(yùn)動模態(tài)和控制面偏轉(zhuǎn)。

2.2 廣義非定常氣動力模型

廣義非定常氣動力建模是建立飛行器氣動伺服彈性系統(tǒng)狀態(tài)空間方程的重要步驟。目前工程上常采用有理函數(shù)擬合方法進(jìn)行廣義非定常氣動力建模。主要有以下3種擬合方法[7,8]：

a）最小二乘法（Least Squares，LS）；

b）最小狀態(tài)法（Minimum State，MS）；

c）修正的矩陣 Padé法（Modified Matrix Padé，MMP）。

LS法氣動力擬合計(jì)算量小、精度最高，但氣動增階最大；MMP方法精度介于LS法與MS法之間，而氣動增階也較大；MS法計(jì)算精度略有降低，但仍適合工程應(yīng)用，其所導(dǎo)致的狀態(tài)空間氣動彈性模型階數(shù)最小。

采用LS法進(jìn)行非定常氣動力擬合，可以通過減少擬合模態(tài)數(shù)及滯后根個數(shù)降低被控對象階數(shù)。氣動系數(shù)矩陣Q的LS法有理函數(shù)近似式Qap為[6]

式中 rj為實(shí)常數(shù)，稱作滯后根。

式中 b為參考長度；ω為振動圓頻率；σ為運(yùn)動衰減率；V為飛行速度。

于是可將彈性飛行器的線性小擾動方程轉(zhuǎn)化為狀態(tài)空間形式：

式中 AS為系統(tǒng)矩陣；xS為狀態(tài)向量；BS為控制矩陣；uS為輸入向量。

2.3 傳感器模型

飛行器常采用加速度計(jì)和角速率陀螺來感受機(jī)體運(yùn)動信息，這些傳感器測得的信號是飛行器相對于慣性參考系的速度和角速度在儀器敏感軸上的分量[3]，這些運(yùn)動信號受儀器位置和彈性振動的影響。輸出方程可以表示為

式中 CS為輸出矩陣；DS為直接傳遞矩陣。

2.4 伺服作動器模型

飛行器采用舵機(jī)來驅(qū)動控制面偏轉(zhuǎn)，常用的有液壓舵機(jī)、電動舵機(jī)和氣動舵機(jī)?？紤]常用的電液舵機(jī)，其舵機(jī)傳遞函數(shù)一般可表示為[3,9]

式中ec為輸入控制指令；δ為控制面偏轉(zhuǎn)響應(yīng)。

將此傳遞函數(shù)進(jìn)行變換，可得到以 ec為控制向量的舵機(jī)狀態(tài)空間方程。

2.5 氣動伺服彈性閉環(huán)系統(tǒng)模型

將彈性機(jī)體環(huán)節(jié)、非定常氣動力環(huán)節(jié)、傳感器環(huán)節(jié)、舵機(jī)環(huán)節(jié)等聯(lián)接起來，可得飛行器ASE廣義受控對象模型：

式中AP為廣義受控對象系統(tǒng)矩陣；BP為控制矩陣；CP為輸出矩陣； xP為廣義受控對象狀態(tài)向量；u為作動器控制指令；y為傳感器輸出的運(yùn)動信息。

這個模型可以用于氣動伺服彈性系統(tǒng)的穩(wěn)定性分析和主動控制律設(shè)計(jì)。

飛行器控制系統(tǒng)的 k階控制律可用狀態(tài)空間形式表達(dá)為

式中 AK為系統(tǒng)矩陣；BK為控制矩陣；CK為輸出矩陣。

將廣義受控對象與飛行控制系統(tǒng)進(jìn)行聯(lián)立，得到飛行器氣動伺服彈性閉環(huán)系統(tǒng)方程：

式中r為參考輸入。

至此，可采用該模型進(jìn)行飛行器ASE系統(tǒng)穩(wěn)定性分析。

3 跨大氣層飛行器氣動伺服彈性穩(wěn)定性分析

3.1 氣動伺服彈性魯棒穩(wěn)定性

對于完全不受攝動的線性系統(tǒng)，其動態(tài)特性由狀態(tài)空間方程描述，穩(wěn)定性由狀態(tài)矩陣的特征值來判斷。但在實(shí)際當(dāng)中，不確定性是客觀存在的。由于不確定性攝動的影響，真實(shí)系統(tǒng)工作時(shí)有可能偏離所預(yù)料的性能，甚至失去穩(wěn)定性。魯棒性（Robustness）即為動力學(xué)系統(tǒng)在一定的不確定性攝動影響下，仍能穩(wěn)定且正常工作的能力。對于跨大氣層飛行器，由于飛行過程中氣動伺服彈性系統(tǒng)模型參數(shù)攝動和外界復(fù)雜擾動因素的影響，氣動伺服彈性魯棒穩(wěn)定性比傳統(tǒng)的穩(wěn)定性更具有現(xiàn)實(shí)意義。

3.2 Nyquist方法

對于如圖1所示ASE閉環(huán)系統(tǒng)，其開環(huán)傳遞函數(shù)可以由各個環(huán)節(jié)的傳遞函數(shù)得到：

在關(guān)心的頻率范圍內(nèi)計(jì)算系統(tǒng)開環(huán)傳遞函數(shù)，在復(fù)平面內(nèi)作圖，可得系統(tǒng)開環(huán)傳遞函數(shù)的乃氏圖，如圖 2所示。由此可以得到閉環(huán)系統(tǒng)的幅值裕度和相位裕度。相位裕度即開環(huán)乃氏曲線上模值等于1的矢量與負(fù)實(shí)軸的夾角，在伯德圖上相當(dāng)于處的相位與π-的角差；幅值裕度即開環(huán)乃氏曲線與負(fù)實(shí)軸相交點(diǎn)模值的倒數(shù)，在伯德圖上，即相當(dāng)于相位為 π-時(shí)，幅值的絕對值，即：

圖2 傳統(tǒng)的穩(wěn)定裕度示意

上述穩(wěn)定裕度概念僅在一定程度上反映系統(tǒng)的魯棒穩(wěn)定性，不能考慮系統(tǒng)是否能承受某類具體的參數(shù)攝動。例如，對于如圖3所示的特定情況[10]，傳統(tǒng)的穩(wěn)定裕度表征方法失去了其有效性。

圖3 傳統(tǒng)穩(wěn)定裕度概念的失效

對于多輸入、多輸出（Multiple Input Multiple Output，MIMO）系統(tǒng)，相應(yīng)的也有廣義Nyquist穩(wěn)定判據(jù)，但該方法可能給出過于樂觀的結(jié)論。在線性系統(tǒng)范圍內(nèi)比較有效的分析工具是最小奇異值方法和結(jié)構(gòu)奇異值μ方法，它們能較好地解決線性系統(tǒng)的魯棒性問題。其中最小奇異值方法只能處理非結(jié)構(gòu)化的魯棒穩(wěn)定性問題，對于參數(shù)不確定性的魯棒穩(wěn)定性分析的問題，比較有效的方法是結(jié)構(gòu)奇異值μ方法。

3.3 最小奇異值方法

對于MIMO系統(tǒng)，以最小奇異值原理為基礎(chǔ)的魯棒穩(wěn)定性分析方法為MIMO系統(tǒng)提供了一種可用的魯棒穩(wěn)定性分析工具。

典型的受乘法攝動的閉環(huán)控制系統(tǒng)如圖4所示。

圖4 反饋系統(tǒng)乘法攝動

攝動矩陣L可表示為

L為氣動伺服彈性系統(tǒng)攝動項(xiàng)的增益kj和相位φj組成的對角矩陣，每一項(xiàng)表示對應(yīng)輸入分量的攝動量。對于標(biāo)稱系統(tǒng)，L為單位陣。

標(biāo)稱系統(tǒng)的回差矩陣為

基于奇異值理論的魯棒穩(wěn)定性判據(jù)為若標(biāo)稱系統(tǒng)穩(wěn)定，則受擾系統(tǒng)穩(wěn)定的充分條件是在全頻域范圍內(nèi)，有以下不等式成立[11]：

實(shí)際計(jì)算中，為滿足系統(tǒng)的裕度需要，一般要求閉環(huán)系統(tǒng)回差矩陣的最小奇異值不小于 0.6，即≥0.6。

3.4 結(jié)構(gòu)奇異值μ方法

基于最小奇異值理論的魯棒穩(wěn)定性分析方法通常給出偏保守的結(jié)果，而結(jié)構(gòu)奇異值μ方法考慮了不確定性攝動集合的塊對角結(jié)構(gòu)，克服了分析中的保守性，能更加精確地對具有結(jié)構(gòu)和非結(jié)構(gòu)不確定性攝動的系統(tǒng)進(jìn)行魯棒穩(wěn)定性分析。

對于線性控制系統(tǒng)，分析其魯棒穩(wěn)定性時(shí)可按照在標(biāo)稱系統(tǒng)中加入攝動影響因子的方法來考慮，如圖5所示。

圖5 魯棒穩(wěn)定性分析

M(s)和Δ(s)均是穩(wěn)定的。在某一頻率處，廣義系統(tǒng)算子M∈Cn×n，不確定性算子Δ是由重復(fù)標(biāo)量塊和滿塊矩陣組成的復(fù)數(shù)塊對角的矩陣，其集合記為

復(fù)數(shù)矩陣M∈Cn×n關(guān)于不確定性Δ的結(jié)構(gòu)奇異值μ(M)為[12]

結(jié)構(gòu)奇異值μ與系統(tǒng)魯棒穩(wěn)定性關(guān)系為：設(shè)α＞0，對于所有滿足的Δ(s)，如圖5所示的閉環(huán)系統(tǒng)是魯棒穩(wěn)定的充分必要條件是[13]：

結(jié)構(gòu)奇異值μ分析方法就是將分析的對象轉(zhuǎn)化為圖 5所示的標(biāo)準(zhǔn)模型，計(jì)算在不確定性攝動矩陣Δ作用下的廣義系統(tǒng)M的結(jié)構(gòu)奇異值 ()μM ，據(jù)此判斷系統(tǒng)的魯棒穩(wěn)定性。

4 算例分析

以某飛行器為例，建立全機(jī)結(jié)構(gòu)動力學(xué)模型及非定常氣動力計(jì)算模型，分別開展結(jié)構(gòu)動特性計(jì)算和非定常氣動力計(jì)算。典型結(jié)構(gòu)模態(tài)及模態(tài)非定常壓力分布如圖6、圖7所示。

圖6 典型結(jié)構(gòu)模態(tài)

圖7 典型模態(tài)非定常壓力分布

首先建立飛行器氣動伺服彈性系統(tǒng)標(biāo)稱模型，在此基礎(chǔ)上，考慮以下不確定因素影響：

a）燃料消耗帶來的全機(jī)總質(zhì)量、轉(zhuǎn)動慣量的變化；

b）氣動加熱效應(yīng)導(dǎo)致的材料彈性模量的下降，以全機(jī)廣義剛度攝動表示；

c）舵機(jī)擬合參數(shù)（截止頻率ωn、阻尼比ξ）的誤差及分散度；

d）控制系統(tǒng)增益K1的攝動。

考慮如表1所示的兩種參數(shù)攝動情況，開展飛行器控制系統(tǒng)三通道（俯仰、滾轉(zhuǎn)以及偏航通道）交聯(lián)耦合的魯棒氣動伺服彈性穩(wěn)定性分析。

表1 參數(shù)攝動情況

分析結(jié)果μ值曲線如圖8、圖9所示。

圖8 情況1的μ值曲線

圖9 情況2的μ值曲線

根據(jù)μ方法魯棒穩(wěn)定性判據(jù)，在第1種參數(shù)攝動情況下，系統(tǒng)最大μ值小于1，系統(tǒng)魯棒穩(wěn)定；在第2種參數(shù)攝動情況下，系統(tǒng)最大μ值幾乎等于1，系統(tǒng)魯棒穩(wěn)定性較差，極易發(fā)生氣動伺服彈性失穩(wěn)，在設(shè)計(jì)中應(yīng)引起重視。

5 結(jié)束語

跨大氣層飛行器氣動伺服彈性系統(tǒng)模型建模涉及多個學(xué)科，包含多個子系統(tǒng)，建模復(fù)雜。此外，飛行環(huán)境變化大，飛行速度跨越亞、跨、超聲速范圍；氣動熱效應(yīng)顯著；由于燃料消耗的原因，系統(tǒng)的質(zhì)量特性變化較大。以上因素不可避免地給氣動伺服彈性系統(tǒng)建模帶來系統(tǒng)誤差和外界干擾。為保證安全飛行，有必要進(jìn)行氣動伺服彈性魯棒穩(wěn)定性分析。

對于線性系統(tǒng)，氣動伺服彈性穩(wěn)定性分析方法有Nyquist方法、最小奇異值法以及結(jié)構(gòu)奇異值μ方法等。研究表明，相比其他方法，μ方法適用于控制通道間存在交聯(lián)耦合的飛行器進(jìn)行氣動伺服彈性穩(wěn)定性分析，且能考慮系統(tǒng)參數(shù)不確定性和多種不確定性攝動的影響，因而非常適合跨大氣層飛行器，是一種很有前景的飛行器氣動伺服彈性魯棒穩(wěn)定性分析方法。

[1]陳桂彬, 鄒叢青, 楊超. 氣動彈性設(shè)計(jì)基礎(chǔ)[M]. 北京: 北京航空航天大學(xué)出版社, 2004.

[2]Richetts R H, Noll T E, Whitlow Jr. An overview of aeroelasticity studies for the national aero-space plane[R]. AIAA 93-1313-CP, 1993.

[3]惠俊鵬, 楊超. 考慮參數(shù)攝動的飛行器氣動伺服彈性魯棒穩(wěn)定性分析研究[J]. 導(dǎo)彈與航天運(yùn)載技術(shù), 2015(6): 1-5.

[4]Spainin C V, Zeiler T A, Gibbons M D. Aeroelastic character of a national aerospace plane demonstrator concept[R]. AIAA 93-1314-CP, 1993.

[5]Changho N, Youdan K, Jeffrey B L. Active aeroelastic wing design for gust load alleviation and flutter suppression[R]. AIAA-97-1265-CP, 1997.

[6]吳志剛, 楊超. 氣動伺服彈性系統(tǒng)不確定性建模與魯棒穩(wěn)定性[J]. 航空學(xué)報(bào), 2003,24(4): 312-316.

[7]宋晨, 楊超, 吳志剛. 3種氣動彈性狀態(tài)空間建模方法的對比[J]. 航空學(xué)報(bào), 2007,28(8): 81-86.

[8]Tiffang S H, Karpel M. Aeroservoelastic modeling and applications using minimum-state approximations of the unsteady aerodynamics[R]. AIAA paper 89-1188, 1989.

[9]Mukhopadhyay V. Flutter suppression control law design and testing for the active flexible wing[J]. Journal of Aircraft, 1995, 32(1): 45-51.

[10]Balas G J, Doyle J C, Glover K. MATLAB μ-analysis and synthesis toolbox user’s guide[M]. USA: The MathWorks Inc, 1998.

[11]Newsom J R, Mukhopadhyay V. A multiloop robust controller design study using singular value gradients[J]. Journal of Guidance, Control and Dynamics, 1985, 8(4): 514-519.

[12]Doyle J C. Structured uncertainty in control systems design[C]. USA: Proceedings of IEEE Conference on Decision and Control, 1985.

[13]Lind R, Brenner M. Analyzing aeroservoelastic stability margins using the μ-method[R]. AIAA-98-1895, 1998.

Research on Methodology of Aeroservo-elastic Stability Analysis for Trans-atmospheric Vehicle

Hui Jun-peng1,2, Yang Chao2
(1. R&D Center, China Academy of Launch Vehicle Technology, Beijing, 100076; 2. Beihang University, Beijing, 100191)

Using wing-body configuration and big lifting surface and control surface, trans-atmospheric vehicle usually has low overall structure stiffness and wide-band multi-way cross-linked flight control system, and must withstand complex flight environment and serve aerodynamic heating during reentry. Factors above may easily induce aeroservo-elastic stability problems. Due to external disturbance, vehicle model perturbation and coupling of control channels, robustic aeroservo-elastic stability analysis must be applied to trans-atmospheric vehicle. In this paper, motion equations of flexible flight vehicle are deduced by using Lagrangian method considering coupling between rigid-body modes and elastic vibration modes. State-space model of closed-loop aeroservo-elastic system is established based on rational function approximation of unsteady aerodynamics and modeling of servosystem and flight control system. Furthermore, nyquist criteria, the minimal singular value theory and the structured singular value analysis are studied for aeroservo-elastic stability analysis, and some conclusions are drawn consequently.

Trans-atmospheric vehicle; Aeroservo-elasticity; Stability; Structured singular value; μ-method

V411

A

1004-7182(2017)01-0042-06

10.7654/j.issn.1004-7182.20170111

2016-07-20；

2016-10-20；數(shù)字出版時(shí)間：2017-01-06；數(shù)字出版網(wǎng)址：www.cnki.net

惠俊鵬（1981-），男，高級工程師，主要研究方向?yàn)轱w行器總體設(shè)計(jì)