聶雪媛，劉中玉，楊國(guó)偉

（中國(guó)科學(xué)院 力學(xué)研究所 流固耦合系統(tǒng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，北京 100190）

對(duì)飛行器氣動(dòng)彈性不穩(wěn)定性預(yù)測(cè)在設(shè)計(jì)中非常重要。隨計(jì)算流體力學(xué)（Computational Fluid Dynamics，CFD）、計(jì)算結(jié)構(gòu)力學(xué)（Computational Structure Dynamics，CSD）及計(jì)算機(jī)性能的發(fā)展，直接耦合CFD／CSD的數(shù)值模擬方法成為非線性氣動(dòng)彈性分析中可信度最高方法［1］。然而在氣動(dòng)彈性研究中，非定常氣動(dòng)力計(jì)算占據(jù)大部分計(jì)算時(shí)間，使基于CFD的氣動(dòng)彈性模擬難以在實(shí)際工程應(yīng)用。如何既能保證計(jì)算精度又能提高計(jì)算效率成為氣動(dòng)彈性的研究熱點(diǎn)。非定常流場(chǎng)模型降階技術(shù) （ReducedOrder Models，ROMs）的出現(xiàn)為CFD／CSD數(shù)值模擬方法用于工程設(shè)計(jì)成為可能。降階模型由全流場(chǎng)全階CFD模型近似投影獲得低階模型，該模型能捕捉非線性流場(chǎng)動(dòng)力學(xué)特征，計(jì)算效率遠(yuǎn)高于CFD。非定常氣動(dòng)力ROMs代替CFD與結(jié)構(gòu)模型耦合，能進(jìn)行氣動(dòng)彈性分析或氣動(dòng)伺服彈性分析及設(shè)計(jì)［2－3］。氣動(dòng)力降階模型基于動(dòng)態(tài)線性假設(shè)［4］，在給定非線性定常流場(chǎng)基礎(chǔ)上利用小擾動(dòng)理論推導(dǎo)獲得到線性模型。在基于系統(tǒng)辨識(shí)的CFD降階模型方法中，Silva等［5］提出基于脈沖響應(yīng)的Volterra級(jí)數(shù)模型降階方法，并與CFL3D結(jié)合進(jìn)行氣動(dòng)彈性顫振分析。但此方法需獲得系統(tǒng)全部Markov參數(shù)，而所用辨識(shí)數(shù)據(jù)因通過“一次激勵(lì)結(jié)構(gòu)一個(gè)模態(tài)”獲得，會(huì)造成結(jié)構(gòu)模態(tài)個(gè)數(shù)較多時(shí)所需脈沖響應(yīng)個(gè)數(shù)呈冪函數(shù)增長(zhǎng)。為此，文獻(xiàn)［6－7］開發(fā)出同時(shí)激勵(lì)結(jié)構(gòu)多個(gè)模態(tài)方法提高辨識(shí)效率。認(rèn)為輸入多個(gè)激勵(lì)信號(hào)間須存在相位差；但未給出如何確定相位差，需多次調(diào)試。Gupta等［8－10］提出采用飛行試驗(yàn)中易實(shí)現(xiàn)的速度為3211信號(hào)作為輸入進(jìn)行系統(tǒng)辨識(shí)。利用自回歸滑動(dòng)平均模型（Autoregressive Moving Average，ARMA）作為替代CFD的降階模型。該方法亦需對(duì)各模態(tài)同時(shí)輸入的3211信號(hào)進(jìn)行反復(fù)試驗(yàn)、調(diào)整頻率，以期輸出響應(yīng)能含希望激發(fā)的頻帶。為進(jìn)一步提高效率，Raveh［11］用多個(gè)白噪聲信號(hào)模擬結(jié)構(gòu)多模態(tài)同步輸入，辨識(shí)三種多輸入、輸出（Multiple Input Multiple Output，MIMO）模型，即頻域模型、ARMA模型及狀態(tài)空間模型，分析AGARD445.6機(jī)翼在跨聲速區(qū)域的顫振邊界。該方法雖可避免對(duì)輸入信號(hào)頻譜范圍要求，但仍未涉及多個(gè)輸入同步信號(hào)間相位差的確定。徐敏等［12－13］以一次激勵(lì)一個(gè)模態(tài)方式的階躍響應(yīng)信號(hào)為輸入，基于Volterra級(jí)數(shù)方法建立非定常氣動(dòng)力降階模型，并用于二維機(jī)翼氣動(dòng)彈性主動(dòng)控制及氣動(dòng)彈性分析。張偉偉等［14］用3211信號(hào)建立基于CFD的 ARMA模型用于 Isogai機(jī)翼及AGARD445.6機(jī)翼氣動(dòng)彈性分析，并認(rèn)為輸入信號(hào)需反復(fù)調(diào)試。劉曉燕等［15］采用基于觀測(cè)器的Volterra核辨識(shí)方法，通過每個(gè)模態(tài)單獨(dú)激勵(lì)氣動(dòng)力系統(tǒng)方式，在每個(gè)模態(tài)輸入含結(jié)構(gòu)感興趣頻率正弦信號(hào)建立氣動(dòng)力降階模型用于氣動(dòng)彈性非線性氣動(dòng)力預(yù)測(cè)及顫振抑制。該方法本質(zhì)仍屬基于Volterra級(jí)數(shù)的ROMs。

為方便輸入信號(hào)選擇、避免同步輸入多信號(hào)時(shí)反復(fù)測(cè)試方可確定相位差、提高辨識(shí)計(jì)算效率，本文以隨機(jī)白噪聲信號(hào)為系統(tǒng)輸入，用一次激勵(lì)一個(gè)模態(tài)方式同時(shí)運(yùn)行多個(gè)CFD數(shù)值模擬程序計(jì)算系統(tǒng)響應(yīng)。利用所得多組訓(xùn)練數(shù)據(jù)辨識(shí)出多個(gè)基于ARMA的單輸入多輸出（Single Input Multiple Output，SIMO）部分氣動(dòng)力降階模型。因非定常氣動(dòng)力建模建立于動(dòng)態(tài)線性假設(shè)理論基礎(chǔ)上，故利用線性系統(tǒng)疊加原理獲得基于CFD的多輸入、輸出非定常氣動(dòng)力降階模型，并將其與結(jié)構(gòu)模型耦合用于AGARD445.6機(jī)翼顫振邊界預(yù)測(cè)。通過將預(yù)測(cè)結(jié)果與風(fēng)洞試驗(yàn)結(jié)果比較表明，本文所提建模方法準(zhǔn)確、有效。

1 基于ARMA的氣動(dòng)力降階模型

顫振邊界預(yù)測(cè)為氣動(dòng)彈性分析中的重要計(jì)算及研究彈性體在小擾動(dòng)激勵(lì)下系統(tǒng)響應(yīng)，符合動(dòng)態(tài)線性模型分類［16］，因此可用線性模型描述非定常氣動(dòng)力，進(jìn)行氣動(dòng)彈性顫振分析。

本文所用氣動(dòng)力降階模型進(jìn)行氣動(dòng)彈性顫振邊界分析過程可描述為將給定的結(jié)構(gòu)模態(tài)位移一次一個(gè)模態(tài)輸入到非定常流場(chǎng)求解器，利用N－S方程求出各階模態(tài)對(duì)應(yīng)的廣義氣動(dòng)力，設(shè)關(guān)心結(jié)構(gòu)的前n個(gè)模態(tài)，則會(huì)獲得n組“單輸入－n輸出”（記為1－n）的模型訓(xùn)練數(shù)據(jù)，用系統(tǒng)辨識(shí)技術(shù)可得n個(gè)“1－n”的線性模型，利用線性系統(tǒng)疊加原理，最終獲得“n個(gè)輸入－n個(gè)輸出”（記為n－n）的氣動(dòng)力降階模型，見圖1下方虛線框。用所得降階模型替代CFD與結(jié)構(gòu)模型耦合，在結(jié)構(gòu)振幅小擾動(dòng)下通過改變動(dòng)壓，可快速獲得結(jié)構(gòu)位移響應(yīng)，從而迅速找到某馬赫數(shù)下的顫振臨界參數(shù)，見圖1上方虛線框。此過程所耗時(shí)間與CFD／CSD耦合計(jì)算相比可忽略。

圖1 基于氣動(dòng)力辨識(shí)模型的氣動(dòng)彈性分析Fig.1 Aeroelastic analysis based on aerodynamic force identification model

本文所用氣動(dòng)力降階模型為自回歸滑動(dòng)平均ARMA模型，即時(shí)間域離散差分方程，形式為

式中：ya為n個(gè)廣義氣動(dòng)力輸出組成的列向量；u為單模態(tài)位移輸入信號(hào)標(biāo)量；na，nb為輸出、輸入延遲階數(shù)；k為離散時(shí)間變量；A，B為欲辨識(shí)系統(tǒng)模型參數(shù)。

設(shè)結(jié)構(gòu)模態(tài)數(shù)為n（以下n均為結(jié)構(gòu)前n階模態(tài)），則式（1）所描述的為圖1下方虛線框中每個(gè)“單模態(tài)位移輸入－n個(gè)模態(tài)廣義氣動(dòng)力輸出”的1－n ARMA模型。據(jù)文獻(xiàn)［14］取4～8階延遲模型即可精確描述非定常效應(yīng)，因此利用MATLAB的arx函數(shù)可較快辨識(shí)最好的表征訓(xùn)練數(shù)據(jù)（式（1））最優(yōu)的1－n ARMA模型。獲得1－n ARMA模型后用文獻(xiàn)［17］方法在Matlab中編程實(shí)現(xiàn)將離散差分方程轉(zhuǎn)換為離散狀態(tài)空間方程形式。第i階模態(tài)激勵(lì)下狀態(tài)空間方程為

上述系數(shù)矩陣均為分塊矩陣。式（2）為第i階模態(tài)激勵(lì)下所得“單輸入－n輸出”的SIMO系統(tǒng)狀態(tài)空間方程。為便于氣動(dòng)彈性分析，需獲得“n個(gè)模態(tài)輸入激勵(lì)下n個(gè)模態(tài)廣義氣動(dòng)力輸出”的多輸入多輸出MIMO系統(tǒng)狀態(tài)空間模型。因此對(duì)n個(gè) “單輸入－n輸出”狀態(tài)空間方程線性疊加，獲得所需MIMO狀態(tài)空間方程為

各參數(shù)中上標(biāo)為模態(tài)階次。

2 激勵(lì)信號(hào)

模型辨識(shí)主要為激勵(lì)信號(hào)選取，要求所用激勵(lì)信號(hào)能激起被辨識(shí)系統(tǒng)所有頻段信號(hào)。隨機(jī)信號(hào)為含頻帶范圍最寬的信號(hào)，理論上用其作為激勵(lì)信號(hào)能獲得系統(tǒng)所有頻段范圍內(nèi)響應(yīng)。為實(shí)際氣動(dòng)彈性工程應(yīng)用，僅研究某特定頻率附近模態(tài)運(yùn)動(dòng)所致氣動(dòng)力，對(duì)隨機(jī)信號(hào)進(jìn)行帶通濾波處理。本文主要研究AGARD445.6機(jī)翼前四階模態(tài) 9.60 Hz、38.16 Hz、48.35 Hz、91.54 Hz引起的廣義氣動(dòng)力。對(duì)輸入信號(hào)進(jìn)行濾波，保證激勵(lì)信號(hào)包含該頻段即可。

輸入信號(hào)類型確定后，數(shù)據(jù)采樣頻率及長(zhǎng)度將決定辨識(shí)精度。由于該激勵(lì)信號(hào)用于CFD計(jì)算，故其采樣時(shí)間與CFD計(jì)算時(shí)間步相同，本文為無量綱時(shí)間0.05，對(duì)應(yīng)有量綱頻率 fs為 11 523.6 Hz（馬赫數(shù)為0.96時(shí)），其 Nyquist頻率為5 761 Hz，遠(yuǎn)高于所需最高頻率91.54 Hz。對(duì)白噪聲激勵(lì)信號(hào)采用帶通方式濾波，保留頻段范圍為0.000 1fs～0.1fs。采樣長(zhǎng)度據(jù)頻率分辨率定義需在CFD上運(yùn)行約2 500個(gè)時(shí)間步，即可辨識(shí)系統(tǒng)最低頻率9.60 Hz。因此可獲得精度較高模型。雖運(yùn)行2 500個(gè)時(shí)間步需耗費(fèi)一定計(jì)算時(shí)間，但與頻繁修改參數(shù)進(jìn)行CFD計(jì)算及氣動(dòng)彈性分析的傳統(tǒng)方法相比，計(jì)算效率改善仍相當(dāng)可觀。

圖2 機(jī)翼近壁面網(wǎng)格圖Fig.2 Mesh grid of near the wing

3 算例與分析

為驗(yàn)證該辨識(shí)方法的有效性及準(zhǔn)確性，本文對(duì)AGARD445.6進(jìn)行模型辨識(shí)，利用辨識(shí)所得模型進(jìn)行顫振邊界、顫振速度預(yù)測(cè)，并與風(fēng)洞試驗(yàn)結(jié)果比較。計(jì)算網(wǎng)格單元數(shù)約40萬。近壁面網(wǎng)格見圖2。CFD網(wǎng)格計(jì)算中首層網(wǎng)格厚度約為根弦長(zhǎng)的8×10－5，結(jié)構(gòu)振型歸一化見文獻(xiàn)［18］。采用非定常 N－S方程計(jì)算機(jī)翼彈性運(yùn)動(dòng)時(shí)非定常流場(chǎng)，空間采用基于結(jié)構(gòu)網(wǎng)格技術(shù)的有限體積法，選HLLE格式進(jìn)行空間離散，用代數(shù)無限插值法實(shí)現(xiàn)網(wǎng)格運(yùn)動(dòng)。用雙時(shí)間步進(jìn)行時(shí)間離散，擬時(shí)間步為隱式格式LUSGS方法。施加模態(tài)位移激勵(lì)時(shí)需保證產(chǎn)生的結(jié)構(gòu)實(shí)際位移、速度滿足線性擾動(dòng)假設(shè)。

本文采用高斯白噪聲信號(hào)，在輸入CFD前先進(jìn)行濾波處理。所選帶通范圍 0.000 1fs～0.1fs，其中 fs為CFD計(jì)算用采樣頻率。對(duì)應(yīng)有量綱帶通約1.152～1 152 Hz。理論上為辨識(shí)出最低頻率，需至少2 500個(gè)時(shí)間步的CFD計(jì)算結(jié)果作為模型訓(xùn)練數(shù)據(jù)，但實(shí)際應(yīng)用中的ARMA模型，用200個(gè)時(shí)間步足夠，多余數(shù)據(jù)可用于校驗(yàn)辨識(shí)的模型。由于針對(duì)結(jié)構(gòu)前四階模態(tài)，因此只要辨識(shí)四個(gè)“單輸入－4輸出”的ARMA模型。模型辨識(shí)結(jié)果階次na＝4，nb＝5。

圖3為馬赫數(shù)0.96、輸入第一階模態(tài)位移為過濾高斯白噪聲信號(hào)時(shí)辨識(shí)模型與CFD計(jì)算所得氣動(dòng)力見圖3，其中虛線為模型計(jì)算結(jié)果，實(shí)線為CFD計(jì)算結(jié)果。兩種方法計(jì)算結(jié)果吻合。

圖4為馬赫數(shù)0.96、給定第二階模態(tài)位移50 Hz正弦信號(hào)時(shí)辨識(shí)所得ARMA模型計(jì)算的氣動(dòng)力與CFD直接計(jì)算的氣動(dòng)力結(jié)果比較見圖4，其中虛線為模型計(jì)算結(jié)果，實(shí)線為CFD計(jì)算結(jié)果，可見二者吻合較好。

圖3 第一階模態(tài)激勵(lì)下辨識(shí)模型與CFD計(jì)算結(jié)果對(duì)比Fig.3 Comparison of ARMA model simulated responses and CFD response to excitation of the 1st mode

圖4 第二階模態(tài)激勵(lì)下辨識(shí)模型與CFD計(jì)算結(jié)果比較Fig.4 Comparison of ARMA model simulated response and CFD response to sinusoidal excitation of the 2nd mode

圖5 輸入3211信號(hào)Fig.5 3211 Excitation of the first four modes

三、四階模態(tài)激勵(lì)的模型準(zhǔn)確性檢驗(yàn)與一、二階類似。在檢驗(yàn)辨識(shí)出的四個(gè)“單輸入－4輸出”ARMA模型準(zhǔn)確性基礎(chǔ)上對(duì)其分別進(jìn)行狀態(tài)空間轉(zhuǎn)換獲得（式（2））方程再進(jìn)行線性疊加，最終獲得“4輸入－4輸出”氣動(dòng)力狀態(tài)空間模型。對(duì)該空間模型精確度采用3211信號(hào)進(jìn)行檢驗(yàn)。輸入具有一定相位差的4個(gè)3211信號(hào)見圖5。在圖5信號(hào)輸入下“4輸入－4輸出”狀態(tài)空間模型計(jì)算的廣義氣動(dòng)力與CFD直接計(jì)算結(jié)果比較見圖6，其中實(shí)線為CFD計(jì)算結(jié)果，方塊為“4輸入－4輸出”狀態(tài)空間模型計(jì)算結(jié)果。由圖6看出，MIMO系統(tǒng)狀態(tài)空間模型計(jì)算結(jié)果與CFD直接計(jì)算結(jié)果吻合較好為評(píng)價(jià)辨識(shí)模型計(jì)算精度，本文用常用評(píng)價(jià)模型準(zhǔn)確度參數(shù)即模型擬合度為

圖6 模型輸出與CFD計(jì)算結(jié)果比較Fig.6 Comparison of identification model simulated results and CFD results

式中：為CFD計(jì)算結(jié)果；為辨識(shí)模型計(jì)算結(jié)果；m為數(shù)據(jù)長(zhǎng)度；η＝1為ARMA模型與N－S方程計(jì)算結(jié)果完全一致。

白噪聲信號(hào)、50 Hz正弦信號(hào)及3211信號(hào)作為結(jié)構(gòu)模態(tài)位移輸入時(shí)，ARMA的模型擬合度見表1。由表1看出，所辨識(shí)模型精度較高。

表1 不同輸入信號(hào)激勵(lì)下模型擬合度Tab.1 Comparison of model fitness to different excitation signals

采用辨識(shí)模型預(yù)測(cè)的AGARD445.6機(jī)翼顫振頻率邊界及速度邊界見圖7、圖8。由二圖看出，耦合結(jié)構(gòu)的N－S方程及ROM計(jì)算結(jié)果與風(fēng)洞試驗(yàn)結(jié)果均較吻合，且機(jī)翼顫振速度在跨聲速區(qū)域出現(xiàn)特有的“凹坑”現(xiàn)象。表明本文所用非定常氣動(dòng)力模型辨識(shí)方法能用于跨聲速氣動(dòng)彈性穩(wěn)定性分析。

圖7 AGARD 445.6機(jī)翼無因次顫振頻率邊界Fig.7 Nondimensional flutter frequency boundary of AGARD 445.6 wing

圖8 AGARD 445.6機(jī)翼無因次顫振速度邊界Fig.8 Nondimensional flutter speed boundary of AGARD 445.6 wing

基于系統(tǒng)辨識(shí)模型進(jìn)行氣動(dòng)彈性分析效率較高。如尋找馬赫數(shù)為0.96時(shí)的顫振臨界參數(shù)。采用求解非定常N－S方程方法搜尋一次參數(shù)至少需CFD運(yùn)行2000個(gè)周期，若每個(gè)周期含20次內(nèi)迭代（1次內(nèi)迭代約3 s），運(yùn)行時(shí)間約需2000 min（Intel Xeon 2.4 GHz雙CPU），而改變參數(shù)重新計(jì)算仍費(fèi)時(shí)2 000 min。用非定常氣動(dòng)力建模技術(shù)時(shí)以產(chǎn)生2 500個(gè)辨識(shí)用數(shù)據(jù)而言，需CFD運(yùn)行2 500 min，而參數(shù)辨識(shí)時(shí)間在Matlab下約需10 s。故模型確定后搜尋一次參數(shù)，Matlab僅需1～2 s。

4 結(jié) 論

（1）本文以過濾隨機(jī)白噪聲信號(hào)為激勵(lì)信號(hào)，利用動(dòng)態(tài)線性假設(shè)及線性系統(tǒng)疊加原理將辨識(shí)所得MISO系統(tǒng)ARMA模型線性疊加獲得MIMO系統(tǒng)模型，實(shí)現(xiàn)非定常氣動(dòng)力建模。將辨識(shí)所得模型用白噪聲信號(hào)、正弦信號(hào)、3211信號(hào)等進(jìn)行檢驗(yàn)，并與直接CFD計(jì)算結(jié)果比較，從而驗(yàn)證辨識(shí)模型的準(zhǔn)確性。該方法與現(xiàn)有辨識(shí)方法相比，具有對(duì)激勵(lì)信號(hào)要求少、易于施加、辨識(shí)效率高等優(yōu)點(diǎn)。

（2）將辨識(shí)所得時(shí)域信號(hào)變換為狀態(tài)空間模型用于AGARD445.6機(jī)翼顫振邊界預(yù)測(cè)。氣動(dòng)力降階模型與結(jié)構(gòu)模型耦合計(jì)算結(jié)果與CFD／CSD耦合計(jì)算結(jié)果及風(fēng)洞試結(jié)果一致，表明本文方法不僅效率顯著提高，能用于氣動(dòng)彈性穩(wěn)定性分析，并為氣動(dòng)伺服彈性分析及設(shè)計(jì)提供可能。

［1］陳剛，李躍明.非定常流場(chǎng)降階模型及其應(yīng)用研究進(jìn)展與展望［J］.力學(xué)進(jìn)展，2011，41（6）：686－701.CHEN Gang，LI Yueming.Advances and prospects of the reduced order model for unsteady flow and its application［J］.Advances in Mechanics，2011，41（6）：686－701.

［2］Dowell E H.Eigenmode analysis in unsteady aerodynamics：Reducedorder models［J］.AIAA Journal，1996，34（8）：1578－1583.

［3］Lucia D J，Beran P S，Silva W A.Reducedorder modeling：new approaches for computational physics［J］.Progress in Aerospace Sciences，2004，40（1／2）：51－117.

［4］Dowell E H.A modern course in aeroelasticity，3rd Revised and Enlarged Edition［M］.New York：Klewer Academic Publishers，1995.

［5］Silva W A.Bartels R E.Development of reducedorder models for aeroelastic analysis and flutter prediction using the CFL3Dv6.0 code［J］.Journal of Fluids and Structures，2004，19：729－745.

［6］Silva W A. Simultaneous excitation of multipleinput multipleoutput CFDbased unsteady aerodynamic systems［J］.Journal of Aircraft，2008，45（4）：1267－1274.

［7］Kim T，Hong M，Bhatia K G，et al.Aeroelastic model reduction for affordable computational fluid dynamicsbased flutter analysis［J］.AIAA Journal，2005，43：2487－2495.

［8］Gupta K K，Voelker L S，Bach C，et al.CFDbased aeroelastic analysis of the X43 hypersonic flight vehicle［C］.Proceedings of the 39th Aerospace Sciences Meeting and Exhibit，2001.

［9］Cowan T J，Gupta K K.Accelerating CFDbased aeroelastic predictions using systems identification［C］.36th AIAA Atmoshperic Flight Mechanics Conference and Exhibit，1998：85－93.

［10］Cowan T J，Gupta K K. Development of discretetime aerodynamic model for CFDbased aeroelastic analysis［C］.Proceedings of the 37th Aerospace Sciences Meeting and Exhibit，1999.

［11］Raveh D E. Identification of computationalfluiddynamic based unsteady aerodynamic models for aeroelastic analysis［J］.Journal of Aircraft，2004（41）：620－632.

［12］徐敏，陳剛，陳士櫓.基于非定常氣動(dòng)力低階模型的氣動(dòng)彈性主動(dòng)控制律設(shè)計(jì)［J］.西北工業(yè)大學(xué)學(xué)報(bào)，2004，22（6）：748－752.XU Min，CHEN Gang，CHEN Shilu.Design method of active control law based on reducedorder model of unste ady aerodynamics for aeroelasticity［J］.Journal of Northwestern Polytechnical University，2004，22（6）：748－752.

［13］姚偉剛，徐敏.基于Volterra級(jí)數(shù)降階模型的氣動(dòng)彈性分析［J］.宇航學(xué)報(bào)，2008，29（6）：1711－1716.YAO Weigang，XU Min.Aeroelastic numerical analysis via Volterra series approach［J］.Journal of Astronautics，2008，29（6）：1711－1716.

［14］張偉偉，葉正寅.基于非定常氣動(dòng)力辨識(shí)技術(shù)的氣動(dòng)彈性數(shù)值模擬［J］.航空學(xué)報(bào)，2006，27（4）：579－583.ZHANG Weiwei，YE Zhengyin.Numerical simulation of aeroelasiticity basing on identification technology of unsteady aerodynamic loads［J］.Journal of Astronautics，2006，27（4）：579－583.

［15］劉曉燕.基于非定常氣動(dòng)力降階模型的氣動(dòng)彈性研究［D］.北京：北京航空航天大學(xué)，2011：54－61.

［16］Dowell E H，Crawley E F.A modern course in aeroelasticity［M］.Kluwer Academic，1989.

［17］Cowan T J， Gupta K K. Development of discretetime aerodynamic model for CFDbased aeroelastic analysis［R］.AIAA，1999.

［18］Carson Y E.Agard standard aeroelastic configurations for dynamic responseiwing 445.6［C］.AGARD－R－756，1988.